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ABSTRACT 


As the demand for Unmanned Aerial System (UAS) technology increases, the current 
guidance, navigation, and control (GNC) algorithms will scale poorly to meet the demand 
because currently, significant resources are required to certify flight controllers on an 
individual platform basis. As different airframes are introduced to meet the expanding 
mission requirements, the resources required to sustain the GNC certification demand will 
become a limiting factor in scalability. The feasibility of replacing conventional GNC 
techniques with modern adaptive control theory was conducted on a commercial-off-the- 
shelf (COTS) open-source autopilot. This enabled rapid prototyping and integration of an 
adaptive controller. The adaptive controller architecture was designed to be aircraft non¬ 
specific. This ensures the controller easily integrates into any aircraft, therefore minimizing 
the resource burden of tuning and certification. The adaptive controller tested in this 
research improved performance over the baseline controller and was rapidly integrated on 
multiple various airframes with minimal resources. Improved performance over classical 
feedback was achieved with fast and robust adaptation in multiple regimes of flight. 
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Executive Summary 


The primary objective of this research was to determine if modern control techniques could 
provide engineering cost and schedule savings for DOD/NAVY autonomous systems. A 
waterfall systems engineering technique was utilized to evaluate the use of adaptive control 
on fixed-wing unmanned aircraft. The growing demand for unmanned systems will inherit 
the costs associated with guidance, navigation, and control. With the use of modern control 
techniques, these costs could potentially be reduced—if not eliminated—as well as gaining 
improved performance over the classical methods. 

The field of adaptive control offers techniques for increasing performance and robustness 
in numerous settings and applications. Adaptive control is different from traditional feed¬ 
back in that it provides a mechanism for adjusting the controller’s parameters to reduce 
plant uncertainty. Classical feedback control utilizes parameters, which are specified by the 
engineer to optimize an ideal use case, which often requires extensive tuning and testing. 
Adaptive controllers adjust their control parameters using various intelligent mechanisms 
designed to increase robustness to plant variation or unanticipated disturbances. Adaptive 
control has many applications in the aerospace domain to include control strategies when 
aerodynamic coefficients are unknown or are non-constant, actuator failure, airframe dam¬ 
age, etc. This research evaluates fixed wing UAS controller performance and robustness 
using the Xi adaptive control architecture. 

Successful integration of the Xi adaptive control algorithm was achieved utilizing dis¬ 
cretization techniques on the Pixhawk 1 commercial autopilot. This algorithm was then 
prototyped and tested utilizing MavProxy and X-Plane 10 as a complete software in the loop 
tool-chain. This enabled rapid prototyping and testing of the Xi algorithm in a high fidelity 
environment ensuring successful integration onto prototype aircraft. 

Two different prototype aircraft were constructed to perform flight tests of the Xi adaptive 
controller with respect to the main objective of reducing the engineering demand required 
to achieve successful flight. Multiple flights were conducted across both prototype aircraft 
and achieved a drastic reduction the effort to achieve fully autonomous capabilities that also 
outperformed conventional PID controllers. 


xvn 




Overall, the Xi adaptive controller implemented in this research met the primary objective of 
reducing the engineering demand required of guidance navigation and control architectures. 
This key performance benefit is realizable through the use of modern control techniques, 
which also resulted in increased controller performance as well as improved battle damage 
tolerance as ancillary byproducts of the improved architecture. Utilizing modern adaptive 
control techniques has been shown as a reliable method for improved performance and 
robust control which could result in immediate cost savings to the DOD/NAVY. 
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CHAPTER 1: 
Introduction 


1.1 DOD / NAVY Autonomous Roadmap 

The Department of Defense (DOD) eondueted an analysis of the role of autonomy, whieh 
outlined technology gaps and predicted advancements required to meet the growing perfor¬ 
mance demands [1]. The study amplifies the fact that autonomy is a challenging field and 
that it is arguably in its infancy. The roadmap attempts to guide decision makers in ensuring 
capitalization of under-utilized technology and succinct awareness of technical challenges 
limiting the current state of the art. Figure 1.1 outlines these elements at increasing scope 
of control comprised of various technology portfolios. The language referenced in the 
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Figure 1.1. DOD Autonomy Roadmap. Source [1]. 

roadmap often recommends that machine learning and/or artificial intelligence is needed 
for elements such as “Fault Detection.” On the lowest level, the roadmap annotates the 
use of guidance navigation and control (GNC) as neither needing improvement or current 
technology being underutilized. An in-depth look at the current state of the art with respect 
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to GNC reveals that GNC is still very costly and arguably antiquated. Research of adaptive 
control as applied to GNC offers new strategies offering improved performance and future 
reduced cost. 


1.2 Challenges in Designing Versatile Controllers 

The unmanned aerial system (UAS) has evolved tremendously over the past decade. Minia¬ 
ture autopilots have gotten smaller and cheaper with more sensitive and redundant sensor 
packages largely due to the cellular phone industry accelerating microelectromechanical 
systems (MEMS) technology. The ability to manufacture these autonomous systems at 
fractions of the cost enables the advancement in multiple cooperative UAS applications in¬ 
cluding swarming capability. This ability to mass-produce large quantities of UAS’s poses 
an interesting challenge. Even though the price has gone down and the performance has 
gone up, there still exists a significant amount of man-hours dedicated to sensor calibration 
and autopilot control law configuration and tuning for optimal performance. 

Over the past decade UAS avionics have drastically improved, but the fundamental control 
law algorithms have not changed. The proportional integral derivative (PID) architecture 
found its origin in automatic ship steering applications in 1922 [2]. Conventional control 
law architectures for UASs predominately still use PID controllers. Their architecture offers 
a well understood and predictable behavior for the class of linear systems. Eor this reason, it 
is well suited for the aviation application. The detriment of PID control is that its application 
is mostly constrained by its use on a linear plant and most aerospace applications are non¬ 
linear and time varying. An aircraft’s control authority that increases proportionally to 
dynamic pressure is one example of significant changes in aerodynamic non-linear control 
behavior. In this case, the PID controller’s robustness to changes in velocity and/or density 
altitude is not guaranteed and for most aircraft has to be delicately handled with lookup 
tables produced from hours of flight test for given configurations. 

Conventional controllers (PID) are difficult to tune and achieving an adequately tuned 
controller requires a significant amount of time and resources. These difficulties can arise 
because of many uncertainties with respect to the aircraft as outlined in the following three 
subsections. 
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1.2.1 Unknown Constant Airframe Parameters 

In the case where airframes are assembled by the same manufacturer, all aircraft still 
require a tedious quality assurance check. Physical aspects of the airframes such as center 
of gravity (CG), control surface deflection/calibration/speed, and airframe alignment all can 
drastically vary within the same delivered batch of airframes. Additionally, these miniature 
UASs experience hard landings, crashes, and/or damage in transportation, which all can 
affect aerodynamic handling qualities. 

1.2.2 Unknown Non-Constant Airframe Parameters 

Aircraft with large airspeed envelopes or various configurations (vertical take off and land¬ 
ing (VTOL), flaps, retractable landing gear), exhibit changing dynamics that challenge 
conventional PID robustness. Because these dynamics are changing drastically, the PID 
controllers either underperform or have to be over designed with ad-hoc gain scheduling 
techniques which significantly increases the controller’s already complex tuning process. 

1.2.3 Fault Tolerance 

If an aircraft exhibits airframe faults or battle damage, the desirable outcome results in the 
controller achieving robust performance. Classical PID controllers have limited to no ability 
to cope with these types of failures. Conventional methods require a priori knowledge of 
the failure scenario. Reconfiguration schemes are developed for specific failures, further 
complicating designs. 

1.3 Problem Formulation and Thesis Organization 

This thesis addresses a problem of shortening the time for developing and fielding a new 
autonomous aerial vehicle by utilizing an adaptive controller, which relaxes some of the 
traditional constraints and assures better robustness in the case of airframe configuration 
variation, various uncertainties, and faults. 

To address the formulated problem, this thesis is organized as follows. 

Chapter 2 provides an overview of modern adaptive control. The brief history of adaptive 
control is introduced to further amplify the specific use case for aerospace applications. 
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The adaptive control architecture is compared with conventional feedback controllers (PID) 
to clarify the distinction between the two approaches for stable and robust control. The 
model reference adaptive control (MRAC) architecture is defined in order to articulate the 
difference between the MRAC architecture and the specific modifications utilized in this 
research to formulate the Xi adaptive controller. 

Chapter 3 outlines the proposed adaptive controller architecture used instead of conventional 
less-robust controllers. It addresses why the architecture must take the form of indirect 
adaptive control vice direct adaptive control. The definition of estimated parameters are 
defined with respect to the aerodynamics model derived in Appendix B. The filter section 
of the controller used to bandwidth limit the controller is defined with respect to the specific 
three parameter model chosen for this research. The Xi adaptive controller discretization 
process is reviewed with respect to methods implemented to achieve integration on a 
commercial-off-the-shelf (COTS) autopilot, which is described in Chapter 4. 

Chapter 4 discusses the COTS autopilot and flight stack code chosen for this research. It 
also covers the ground control station (GCS) and software in the loop (SITL) infrastructure 
which was heavily utilized for initial code development, testing, and validation. Chapter 4 
follows with the major steps undertake to prototype two airframes to be used in the flight 
testing of the proposed controller. 

Chapter 5 and 6 describe the results of computer simulations and flight tests, respectively. 
Chapter 6 also outlines the shortcomings of the COTS autopilot architecture and improve¬ 
ments specific to the Xi adaptive controller. 

The thesis concludes with Chapter 7 providing conclusions. 
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CHAPTER 2: 

Overview of Modern Control Techniques 


This chapter is a general review of the history of adaptive control to include its use cases and 
previous pitfalls. A brief overview of algorithm differences between conventional feedback 
control and MRAC architectures is also covered. 


2.1 Adaptive Control History 

Adaptive control saw its early debut in the NASA North American X-15 hypersonic rocket- 
powered X-plane experimental aircraft. The X-15’s performance envelope exceeded Mach 
6.0 and 300,000 feet [3]. Engineers realized early on that the linear controllers performed 
well only at one dynamic pressure, but nowhere near the entire flight envelope. Scheduling 
the controller gains with respect to dynamic pressure (gain scheduling) was one method 
used to help ensure robustness; the method is still widespread in commercial aviation due 
to its robustness but requires a significant effort to explore the entire flight envelope. These 
non-trivial efforts were what encouraged the exploration of the benefits of adaptive control. 

The X-15 program started in 1959 and continued to 1968 flying nearly 200 successful 
flights [4]. It was considered one of NASA’s most successful programs. The benefit of 
adaptive control to the X-15 was that the adaptive controller was supposed to adjust the gain 
parameters online automatically. If the controller was self-tuning, it could potentially offer 
increased performance while reducing complexity. Honeywell implemented the MH-96 
adaptive controller in the X-15-3 as a fly-by-wire controller designed to adaptively adjust 
the damping in pitch and roll with respect to the desired model response. The goal was 
to achieve consistent aircraft response regardless of dynamic pressure and other variables. 
Dydek et. al comments that during test flights of the MH-96 adaptive control, increased 
performance was observed especially in the dynamic phases of re-entry over that of the 
linear fixed gain damping system [5]. These early breakthroughs in adaptive control proved 
the benefits could be viable aerospace solutions. However, on November 15, 1967, there 
was a fatal accident caused by the adaptive controller. NASA documented that the adaptive 
controller created an out of control flight situation resulting in dynamic pressures exceeding 
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the structural limits and subsequent breakup of the airframe at 65,000 feet [4]. 

The turbulent start of adaptive control, as implemented on the X-15 program, was due in 
large part to the early naive understanding of robustness. Contemporary robust adaptive 
control strives to encapsulate these deficiencies of robustness in studies and proofs using 
Lyapunov stability analysis. In addition to the developments of rigorous stability analysis 
tools, a number of unique techniques have also been implemented to increase controller 
robustness. One such technique utilizes dead band limits on the model adaptation process 
to avoid system/measurement noise from causing the un-learning of the states and is called 
“dead-zone modification” as proposed by B.B. Peterson and K.S. Narendra [6]. Lavretsky 
and Wise also reference the “cr and e modifications” which adds damping to the adaptation 
process [7]. The Xi adaptive control algorithm utilizes a technique which seeks to decouple 
the adaptation rate from robustness by “low-pass filtering” the contribution of the fast 
estimator under the premise that estimating the entire frequency spectrum is overly ambitious 
and should be limited to the bandwidth of the actuator-plant combination. Many advances 
have been made in the adaptive control field over the past few decades, and this research sets 
out to evaluate a small subset of these techniques in the unforgiving aerospace environment. 

2.2 Classical Feedback vs Adaptive Control 

Control of a system can be categorized into two required elements; the requirement to 
stabilize the system in the presence of: 

1. disturbances that affect the controlled states and outputs (pitch rate perturbation 
caused by environmental effects) 

2. disturbances that affect the dynamics of the open loop plant (pitch rate effectiveness 
with respect to dynamic pressure) 

Classical feedback control seeks to resolve disturbances that affect the tracking with respect 
to state perturbation and assume that the plant dynamics are constant. This form of control 
is meticulously tuned to achieve the desired overshoot and settling time for example. The 
important assumption that is made by classical feedback controllers is that the underlying 
plant/system dynamics are not changing. For example, the cruise control that maintains 
a vehicle’s speed assumes that the available horsepower of the car is fixed. This is a 
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fairly good assumption as the horsepower with respect to rpm available at sea level and 
5,000 feet for an internal combustion engine is constant enough that a fixed gain feedback 
controller would perform well at maintaining the speed of the vehicle in both environments. 
In the case of an airplane, the dynamic pressure is proportional to velocity squared and 
can drastically change the performance of the aircraft. In this case, the constant system 
performance assumption can cause a fixed gain classical feedback controller to go unstable 
at higher dynamic pressures (higher airspeeds). Conversely, adaptive controllers assume 
that the system performance is unknown and is likely to vary with time. Adaptive control 
seeks to ensure a system’s performance with respect to characteristics, such as damping 
ratio and settling time, which are kept constant regardless of a plant’s dynamics that may be 
unknown and time varying. For both classical and adaptive control, there exists some form 
of error which drives the controller. In the case of classical feedback, the error is calculated 
between the command and the feedback state of the plant. In adaptive control (in general), 
the error is calculated between the outputs of the desired reference model and real plant’s 
measured performance. 

Figure 2.1 outlines the decision making process a controls engineer makes when deciding 
the type of controller needed for a given circumstance. 



Figure 2.1. Determine If Adaptive Control Should Be Used. Adapted from 
[ 8 ]. 
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2.3 Model Reference Adaptive Control 

MRAC establishes the foundation for most of modern, robust adaptive eontrol. Its strueture 
is intuitive in nature and seeks to define a system’s response to a eommand signal with 
a referenee model. Unlike traditional feedbaek where the error signal is generated with 
respeet to state or output error, MRAC’s objeetive is to minimize the error between the 
performanee eharaeteristies of a ehosen referenee model and the plant. The error between 
the model response and the system response generates error for an “adjustment meehanism” 
to learn the unknown model parameters. 

Figure 2.2 illustrates a topology where a traditional feedbaek eontroller is established as 
an inner loop and the “Referenee Model” and “Adjustment Meehanism” is established as 
an outer loop. The outer loop attempts to minimize the error between the referenee model 


reference 

command 



Figure 2.2. Traditional Model Reference Adaptive Control (MRAC) Archi¬ 
tecture. Adapted from [8]. 

output and the plant output. Two primary methods for using this error to learn the system 
parameters are: gradient deseent and Lyapunov stability theory. 

2.3.1 MIT Rule - Gradient Descent 

One of the first approaehes to the design of MRAC eontrollers was implemented at the 
Instrument Labs at MIT (now known as Draper Labs). The gradient deseent based method 
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was called the “MIT Rule” for this reason [8]. This method attempts to learn some unknown 
parameter by descending the gradient of the error between the reference model and the plant 
output. 

Given the simple first-order system G{s): 

G(s) = kdc-^ (2.1) 

5 - 1-1 


where kdc is some unknown feedforward gain. In the case of the MIT rule, kdc is the 
parameter to be learned and is defined as 6. The first step in the MIT rule is to establish a 
cost (or loss) function. One example of a cost function J(6) is: 

J(0) = ( 2 . 2 ) 

e = y - ym (2.3) 


where e is error, y is the plant output, and ym is the model output. 


In order for the cost function to be minimized, the negative gradient of the cost function is 
calculated and used to correct the a priori estimate. This method takes the following form 
where y is the adaptation gain: 


do dJ de 
Tt " 


(2.4) 


The stability of this method is very system dependent and heavily relies on trial and error to 
ensure the adaptation gain (y) is not too high. This usually requires low adaptation rates for 
most systems and may not produce adequate results. It should also be noted that this method 
relies on adequate persistence of excitation [8]. In simple terms, the persistence of excitation 
condition guarantees some functional correspondence between the error and the unknown 
parameters 6. If the error signal is sufficiently rich, then the partial derivative is explicitly 
evaluated by the adaptation process running recursively. Additionally, gradient descents 
are subject to local extreme convergence for non-convex manifolds. Therefore, this method 
does not explicitly guarantee parameter convergence even though correct steady-state gain 
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is achieved for the closed loop system. 

With the differences between classical feedback control and adaptive control being dis¬ 
cussed, the next chapter applies analytical tools to develop a specific adaptive control 
architecture to address the requirements outlined in Chapter 1. 
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CHAPTER 3: 

Engineering of Adaptive Control 


This chapter outlines the pertinent aspects of the Xi adaptive control architecture used for 
this research. The distinction between direct and indirect adaptive control is made in order 
to illustrate the importance of why the indirect architecture is critical for the Xi algorithm. 
This chapter also defines the parameter estimates provided by the Xi algorithm with respect 
to the aerodynamic model simplifications outlined in Appendix B. Finally, the Xi adaptive 
control filter is explained and the methods used for discretizing the algorithm for autopilot 
integration are discussed. It should be noted that the versions of the Xi adaptive controller 
have been patented and proper licensing should be considered if commercial use is desired. 

3.1 £i Adaptive Control 

The Xi adaptive controller is an evolution of the concepts implemented by MRAC. They are 
similar approaches designed to control a system with unknown parameters. The estimates 
of unknown parameters are adjusted to achieve the desired outcome of the error between the 
actual plant (system) and the referenced system model (state predictor) to asymptotically 
approach zero. Adaptive control attempts to estimate the plant’s unknown parameters in 
situ. Parameter estimation is done using either direct or indirect architecture. The indirect 
architecture attempts to estimate the system’s parameters and can be compared to online 
system identification. Alternately, the easier to implement direct architecture estimates the 
controller parameters explicitly. These architectures can be seen in Figures 3.1 and 3.2. 

3.1.1 Reference Model versus Companion Model 

Traditional MRAC controllers often refer to the system objective function as the “Reference 
Model.” In this case, the engineer designs a reference model response, and it is from this 
model response that the error state is calculated directly. Because the Xi adaptive control 
implements the use of a filter in conjunction with a model objective function, the model is 
often referred to as the “Companion Model.” This subtle distinction is necessary because 
the engineer must be aware that the system response will be with respect to the filter plus 
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Figure 3.1. Direct MRAC Architecture 
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Figure 3.2. Indirect MRAC Architecture 


the companion model in series. In other words, the plant will not mimic the companion 
model; it will mimic the companion model plus the filter section. 


3.1.2 Xi Architecture 

The Xi adaptive control algorithm asserts that trying to control the plant uncertainties 
outside of the control actuators’ bandwidth is overly ambitious. The system’s actuator 
bandwidth and the slow dynamics of the plant are most commonly the system’s limiting 
factors, and the estimator’s robustness/stability could be in question if unmodeled high- 
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frequency content exists in the plant. The Xi adaptive control constrains the objective 
function by using a low-pass filter (first or second-order) to limit the frequency response to 
meet robustness specifications. This low-pass filter should be tuned to a frequency response 
commensurate with the actuator’s frequency response. Through inspection the low-pass 
filter placement in the controller topology, it becomes clear that the indirect architecture is 
the only candidate. Figures 3.3 and 3.5 illustrate the placement of the low-pass filter and its 
implication on the closed loop model. 


reference 

command 



Figure 3.3. Direct MRAC Architecture with Low-Pass Filter 

It can be seen in Figure 3.3 that the direct architecture implementation of the low-pass filter 
introduces difficulties. The companion model (state predictor) is defined using aerodynamic 
model and the stability is verified for a Lyapunov candidate function. However, these 
underpinnings designed in the companion model become nonsensical in the implementation 
found in Figure 3.4 because the filter is applied to the plant and not modeled in the companion 
model. This specific filter placement is non-subtractable when the error state is calculated 
because it is not applied in both signal paths (plant and state predictor). 

Conversely, the indirect approach offers an implementation which ensures the low-pass filter 
is applied to both the companion model (state predictor) and the plant ensuring that the 
interaction of the filter is subtractable when calculating the error state. 

It can be seen that the low-pass filter in the direct architecture inherently changes the 
structure of the plant with the cascading of the low-pass filter and plant block diagrams. 
This change mathematically is not mirrored in the companion model (state predictor) and 
therefore is not subtractable. This means that the state predictor is modeling incomplete 
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Figure 3.4. Non-Subtractable Low-Pass Implementation - Direct Architec¬ 
ture 
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Figure 3.5. Indirect MRAC Architecture with Low-Pass Filter 

information and therefore will ineorreetly eharaeterize the system. However, in the indireet 
ease, the strueture of the model is kept intaet, and the low-pass filter is applied to both 
the plant and the state predietor. This ensures that the low-pass filter is subtraetable when 
ealeulating the error state and the model’s strueture is kept intaet. 

Many variations of the Xi adaptive arehiteetures have been derived for various use eases [9]. 
Some of the following forms were studied for viability in the fixed wing UAS use ease with 
respeet to state parameters {6, co, cr) as deseribed in Figure 3.6 and Seetion 3.2: 

• single input single output (SISO) with eonstant but unknown state parameters 
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• SISO with time variant and/or nonlinear unknown bounded state parameters 

• multiple input multiple output (MIMO) with constant but unknown state parameters 

• MIMO with time variant and/or nonlinear unknown bounded state parameters 

MIMO control algorithms would potentially afford the controller more ability to cope with 
system coupling if present. Fixed wing UAS equations of motion, as seen in Equation B.l 1, 
exhibit coupled behavior both in the aerodynamic and inertial coupling. However, MIMO 
was not chosen due to the added complexity required to architect the algorithm into source 
code. Matched uncertainty are terms that can be factored from the control matrix (B) and can 
be time invariant. Unmatched uncertainties are typically modeled as process noise impacting 
the state (i) directly regardless of control. The unmatched uncertainty architecture offers a 
more appealing solution for fixed wing use cases (asymmetric actuator failure, aerodynamic 
coefficients scaled by dynamic pressure), but adds a significant amount of complexity to the 
architecture. In summary, the SISO architecture with matched uncertainty was chosen for 
this research. 

The SISO controller with matched uncertainty was selected to control pitch rate (q) and 
roll rate (p) of the aircraft using two separate controllers. The simplified equations in 
Equation B. 13 assume that there is no aerodynamic or inertial coupling in order to simplify 
the model. This simplification was chosen to make the controller as airframe agnostic 
as possible while enabling very simple applications in code. Utilizing this simplified 
architecture may be at the detriment of controller performance, but as seen by the flight test 
results in Section 5.2, the baseline performance is adequate to achieve performance which 
is better than PID in some cases. In this implementation of the Xi adaptive controller, the 
desired state v to be controlled was an individual body rate (e.g., q, p). 

As seen in Eigure 3.6, the generalized Xi architecture in block diagram form and the 
following elements can be identified 

kg - feed forward input gain 
k - feedback gain 

D{s) - user described filter (second-order low-pass plus integrator) 
rj - JLi controller state 

X - first-order differential equation of state model 
X - state estimate 
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Figure 3.6. Xi Architecture with Matched Uncertainty Block Diagram. 
Source [9]. 


X - state error 
u - eontroller output 
r - referenee input 
Am - Hurwitz matrix 
b - input matrix 

a) - unknown input gain eoeffieient 
0 - unknown eonstant state eoeffieient 
O' - unknown disturbanee estimate 
r - adaptation gain 

Pb - solution to the Lyapunov stability equation 

It should also be noted that the arehiteeture presented in Figure 3.6 ineludes the use of a 
projeetion operator. The estimation dynamics of d), 9, and a are projection based adaptation 
laws. This ensures that the adaptation stays bounded around the feasible region of parameter 
space. The Lyapunov stability proofs for this architecture rely on this method to guarantee 
stability [9]. More discussion on the specific implementation of this operator can be found 
in Appendix E. 

One of the main benefits of using the SISO architecture is in its simplicity so that the 
solution to the Lyapunov stability equation {Pb) utilized in the projection based adaptation 
laws is greatly simplified. 
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In this case, Pb reduces to: 

Pb=:^ (3.1) 

where a)„ is the natural frequency in rad/s for the first-order companion model in discrete 
recursive form assuming DC gain of 1. The proof for this can be found in Appendix F.4. 

3.2 £i Parameter Estimation 

Three model parameters were evaluated in this research. The first model parameter 9 
establishes the baseline architecture and Lyapunov stability proof. Each additional model 
parameter added to the architecture adds complexity to the architecture and its required 
stability proofs. The general progressive build up of the three parameter architecture 
is outlined in the nomenclature for the estimation of a system with unknown constant 
parameters, input/output disturbances, and an unknown system input gain. 

The Xi adaptive control algorithm is primarily used to estimate unknown constant system 
parameters. These system parameters are defined as 6 as seen in the following model: 

x(t) = Ax(t) + b(u(t) + 6^ (t)x(t)) (3.2) 

The second adaptive element is the estimation of the input/output disturbances (cr). This 
additional parameter is implemented in this research as any unmodeled transient dynamics 
such as aerodynamic/inertial coupling. This model takes the form 

x{t) = Ax{t) + b(u(t) + 6^(t)x(t) + cr{t)) (3.3) 

The last adaptive element used for this research was the estimation of unknown system 
input gain (m). Estimating the unknown system input gain offers the controller the ability 
to estimate the actuator effectiveness for changes in dynamic pressure or failed control 
surfaces. This model takes the form 

x{t) = Ax{t) + b((i>(t)u(t) + 9~'^(t)x(t) + cr(t)) (3.4) 
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The model in Equation 3.4 can be paralleled to the aircraft model derived in Equation B.14. 
As a result, the roll rate model takes the similar form 


p = App + bp {cbp6a + 6pp + d-p) (3.5) 

This final model, with the inclusion of all three estimated parameters (a), 6, o'), established 
the architecture tested in this research. The fundamental assumption that the aerodynamic 
and inertial coupling was negligible and set to zero as outlined in Equation B. 13 is agreeably 
a gross assumption which may prove to be inadequate for stable flight. These assumptions 
presume that the estimated input/output disturbance (&) would be adequate to compensate 
for these unmodeled dynamics and were validated by successful flight test in Section 5.2. 

These derivations of the Xi algorithm guarantees that the error between the model and plant 
asymptotically approaches zero, but this does not imply the constraint that the estimated 
parameters are in fact converging to their real values. The algorithm only guarantees that 
the parameters are bounded and therefore it is common to observe the parameters never 
reaching steady-state. The engineer must ensure that the bounds set on the parameters 
are sufficient for the controlled system to not become unstable. In the discrete form, the 
algorithm can have numerical floating point instability if these projection operator bounds 
are too high and not thoroughly tested in simulation. 

3.3 £i Filter - C{s) 

One of the key features of the Xi adaptive controller is that the robustness of the controller is 
decoupled from the adaptation rate. This is handled in the filter section of the Xi architecture 
which is annotated as C(s). To ensure guaranteed stability of the Xi algorithm, C(s) must 
be verified as strictly-proper stable. With the architecture that includes the system input 
gain (co), one cannot simply apply a stand alone filter. The inclusion of co in the architecture 
block diagram in Eigure 3.6 slightly modifies the signal output and takes the following form 

u(s) = -kD(s)(r](s) - kgr(s)) (3.6) 
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where D(s) is the new user defined filter and C(s) now takes the form 


CO e 0.0 = [m/o, 

ojkD(s) 


C(s) = 


(3.7) 


1 + a)kD(s) 


In the ease whieh the user defined funetion D(s) is a simple integrator, C(s) takes the form 


k)(s) = - 
s 

C(s) = 


cok 


(3.8) 


s + (x>k 


The feedbaek gain k should not be assumed to be 1 in this case. As seen in Equation 3.8, the 
C{s) transfer function is a function of k that can have significant influence on the controller 
output. In actual implementation, k was found to be one of the most influential gains in 
the architecture and extremely critical in specifying the transition between robustness and 
performance. The feedback gain k must be set to low enough to prevent high frequencies 
to the actuators, and high enough to ensure the Xi norm stability conditions [9]. 


3.4 Xi Discrete Time Implementation 

Implementing any algorithm on actual autopilot hardware will inevitably force some if not 
all parts of the algorithm to be discretized. Autopilots like the Pixhawk operate at some 
scheduled loop rate for executing the litany of subprograms that measure sensors, calculate 
navigation commands, and much more. In the case of the Pixhawk autopilot, the main loop 
can run up to 400 Hz. At a 400 Hz sampling rate, there is a significant insurance that the 
vehicle’s dynamics bandwidth will be completely defined. However, the ArduPilotMega 
/ Multi-Platform Autopilot (APM) flight stack records all logged parameters also at this 
loop rate and can create log files larger than are reasonably desired. There are a myriad 
of other reasons why the engineer would not want to run at high loop rates, but successful 
flight at the lowest (default of 50 Hz) is desired if adequate performance of the adaptive 
control can be achieved. Failures in early adaptive control were largely in part due to a 
very naive understanding of robustness. Brian Anderson concludes that “it is clear that 
the identification time scale needs to be faster than the plant variation time scale, else 
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identification cannot keep up” [10]. 


The Xi architecture does not require the algorithm to run in sync with the sensor measure¬ 
ments. Ideally, the autopilot main loop rate would remain at its default value (50 Hz), and the 
adaptation section of the Xi controller would run as fast as the CPU would allow. However, 
the APM architecture does not lend itself well to this scheme in its current configuration 
without significant modification to the code base. Therefore, the initial tests presumed that 
the Xi algorithm is running at the same frequency as the main loop of the autopilot. The 
fundamental derivation of the Xi algorithm proves that performance and adaptation gain 
are maximized as the loop time {dt) approaches zero. Increased performance would come 
from higher adaptation loop rates, but the primary focus was to implement the Xi algorithm 
with minimal or no detrimental effects to the current APM codebase. 


3.4.1 Digital Bi-Quad Filter 

The Xi adaptive control algorithm utilizes two specific elements that will require careful 
discretization; the companion model and the low-pass filter. The digital bi-quad filter offers 
a very versatile and straightforward method for accurately implementing the companion 
model and the low-pass filter discretely using its recursive nature. It is a second-order filter 
which uses a finite impulse response (FIR) front end and an infinite impulse response (HR) 
back end requiring four total memory blocks. This topology allows the designer to create 
numerous types of filters (low-pass, high-pass, band-pass) simply by choosing appropriate 
coefficients. If a first-order filter is needed, then the higher order FIR/HR terms can be set 
to zero. Figure 3.7 illustrates this filter’s topology where the FIR structure is the left two 
memory blocks and the HR structure is the right two memory blocks. 

A bilinear Z transform is used to convert the desired S-domain (continuous time domain) 
filter/model into the Z-domain (discrete time domain) to determine the structure of the 
coefficients. 

This derivation for the second-order low-pass model is 


H{s) = 


.2 + ^ + 1 


(3.9) 
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output 


Figure 3.7. Digital Bi-quad Filter Architecture 


where the bi-linear transform eonverts s to z via 




(3.10) 


K is the “pre-warping” faetor whieh aeeounts for the transition of the vertieal s-plane into 
the eireular z-plane as seen in Figure 3.8. 

where coT is 


coT = In ( — 
\Fs 


(3.11) 


K = tan 


tan n— 


(3.12) 
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Figure 3.8. Bi-linear Transform 

Fc is the desired corner frequency of the filter and Fs is the sampling rate (or loop rate of the 
autopilot). This “pre-warping” is critical to ensure that the continuous time cutolf frequency 
desired is correctly established in the discrete implementation. It is the engineer’s discretion 
if pre-warping is required for the appropriate application, but the general guidance is to pre¬ 
warp the Z-domain coefficients if the desired cut-olf frequency is close to Nyquist-Shannon 
sampling theorem frequency {^). It was chosen for this application to always pre-warp the 
coefficients even though the error is small for corner frequencies which are fairly distant 
from Nyquist-Shannon sampling theorem frequency. Continuous calculation of the pre¬ 
warp coefficient was chosen because calculating the tanQ function real time on the CPU 
adds negligible computational strain but offers ease of tuning for the engineer. 

Applying the bi-linear transform to the continuous time second-order low-pass filter results 
in 


H{z) = 


(f)( 


z+l) 


+ 




Q 


+ 1 


(3.13) 


22 












The desired form is 


H{z) = 


bo + b\z ^ + b 2 Z ^ 
ao + a\z~^ + a2Z~^ 


(3.14) 


Reducing Equation 3.13 to match the form in Equation 3.14 results in the following coeffi¬ 
cients 


ao 


ai 


a2 


bo 


= 1 


2(K^ - 1 ) 
+ I + 1 
- I + 1 
+ I + 1 




i^2 + I + 1 


b\ = 2bo 

b2 = bo 


(3.15) 


The bandwidth of the filter Q can be set by the engineer. Eor example, if the pass-band of 
the filter is desired to be flat (Butterworth) then Q can be set equal to Eor this research, 
the following C-H- code segments were used to explicitly calculate the bi-quad low-pass 
filter implementation [11]: 


void DigitalBiquadFilter <T> :: compute_params ( float sample_freq , 
float cutoff_freq , biquad_params &ret) { 
ret . cutoff_freq = cutoff_freq ; 
ret . sample_freq = sample_freq ; 

float fr = sample_freq / cutoff_freq ; 

float K = tanf (M_PI/fr ); //Pre-Warp calculation 

float c = l.Of+2.0f*cosf (M_PI/4.0 f)*K + K*K; 
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re t.bO = K*K/c ; 

ret.bl = 2.0f*ret.b0; 

ret.b2 = ret.bO; 

ret.al = 2.0f*(K*K-1.0f)/c; 

ret.a2 = (1.0 f-2.0 f * cosf (M_PI/4.0 f )*K+K*K) / c ; 

} 

T DigitalBiquadFilter <T> :: apply ( const T &sample , 
const struct biquad_params &params) { 

T delay_element_0 = sample - _delay_element_ 1 * params . al 
- _delay_element_2 * params . a2 ; 

T output = delay_element_0 * params.bO 
+ _delay_element_ 1 * params.bl 
+ _delay_element_2 * params.b2; 

_delay_element_2 = _delay_element_ 1 ; 

_delay_element_ 1 = delay_element_0 ; 

return output; 

} 


This implementation can be employed as the Xi low-pass filter and as the companion model. 
It can be seen in the code segment that K, the pre-warp factor, is explicitly calculated every 
iteration. 

3.4.2 Simplified Bi-quad First-order Model 

In the case of the companion model, a first-order response may be desired. As described 
in equations A.2 and A.4, the discrete first-order model can be derived from a simplified 
Bi-quad as seen in Figure 3.9. It can be seen that the first coefficient of the HR filter is kept 
from this topology. 

The first-order model can be specified by either its time constant (time in seconds to reach 
63% of steady-state) or its -3dB corner frequency. The system takes the form as seen in 
Equation 3.16 when defined by its corner frequency 
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H(s) = (3.16) 

S + COn 

Therefore, the explicit calculation of the Bi-quad coefficients in this case becomes 



bo = I - a\ 


(3.17) 


where a>n is the -3dB corner frequency in radians per second and Fs is the sampling 
frequency measured in Hertz. 

Therefore, the discrete recursive form of the first-order model becomes 


yi+i = aiyi-i + boyi (3.18) 

Another form designed to optimize for speed that is commonly seen in software is 
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float b_0=exp(-f_c/f_s); 
float out+ = (in-out )* b_0 ; 


3.4.3 Euler vs. Trapezoid Rule 

The model estimate, as well as the parameter estimates for the Xi algorithm, are both numer¬ 
ically estimated using discrete integration. The Euler method is a numerical procedure for 
solving ordinary differential equations. The Euler method as applied to discrete integration 
is the fundamental method for recursively integrating a digital signal. The algorithm takes 
the form: 

where h is the uniform step size. 


i/,+1 = yi + hf(ti, yi) 


(3.19) 


The recursive trapezoidal method takes the form 


yi+i = yt + hf(ti, yi) 

h 

yi+\ = yi + yd + /(h+i, ^i+i)] 


(3.20) 


Comparing the accuracy of the two numerical methods for discretely calculating the integral 
of y = can be seen in Eigure 3.10 

As seen in Equation 3.20, the recursive trapezoidal integration method only adds one more 
line of complexity to the algorithm for a significant gain in accuracy and therefore will be 
the chosen method applied for all discrete numerical integration in this research and takes 
the form 


float trap_integration ( float yO, float yl_dot , float dt , float &y0_dot) 
{ 

float yl = yO + (dt/2)*(y0_dot + yl_dot); 
y0_dot = yl_dot; 

return yl; 

} 
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Figure 3.10. Euler vs Trapezoidal Integration Error 


With the Xi adaptive eontrol algorithm defined in diserete form utilizing the arehiteeture 
in Figure 3.6 and the diseretization methods outlined in Seetion 3.4 enabled the simulation 
testing found in Chapter 5. Flight testing was also desired and therefore multiple test aireraft 
were designed and built as seen in Chapter 4 with integrated COTS autopilots. 
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CHAPTER 4: 

Design of Experimental Platform 


This chapter presents the COTS autopilots used in this research as well as outlining the 
Linux toolchain used to conduct |acSITL and GCS operations. This chapter also introduces 
the airframes that were prototyped for testing specific performance characteristics discussed 
in Chapter 5. 

4.1 Pixhawk Autopilot 

The Pixhawk autopilot is a collaborative project among open-source engineers which re¬ 
sulted in a high-performance autopilot which is capable of controlling aircraft, ground 
vehicles, and many others. The primary reason this autopilot was chosen was because of 
the vast amount of support in the developer community. The Pixhawk 1 autopilot hardware 
is effectively obsolete at the time of this writing, but the open-source code base is extremely 
flexible and continues to be ported to new hardware as it becomes available. This has been 
the case for various Raspberry Pi autopilots as well as the Pixhawk 2. The Pixhawk autopilot 
operates using two flight stacks (code base/operating systems); the PX4 flight stack and the 
APM flight stack. This research was implemented on the APM flight stack primarily be¬ 
cause the author’s familiarity with the developer team which offers unparalleled assistance 
to the academic community. The APM codebase also offers a litany of open-source tools 
such as a Linux based GCS, SITL simulator, log analysis tools, and an application program 
interface (API) for the RealFlight 7.5 simulator (high fidelity airframe simulation for small 
aircraft). 

Figure 4.1 illustrates the connection diagram for the Pixhawk 1 autopilot. 

4.1.1 Key Features 

The Pixhawk 1 Autopilot has the following features as found on [12]: 

• 168 MHz / 252 MIPS Cortex-M4F 


29 





* Spektrum DSf/ . . '. ver 

2 Telemetry (radio telemetry) 

3 Telemetry (on-screen display) 

4 USB 

5 SPl (serial peripheral interface) bu> 

6 Power module 

7 Safety switch button 

8 Buzzer 

9 Serial 

10 GPS module 

CAN (controller area network) bus 

12 'C splitter or compass module 

13 Analog to digital converter 6.6 V 

14 Analog to digital converter 3.d V 

15 LED indicator 



reset button 

2 SO card 

3 Flight management reset button 

4 .•-'.;.-.r..uSB port 



*'' : control receiver input 

2S.6us output 

3 Main output- 

4 Auxiliary outpL.‘~ 

Figure 4.1. Pixhawk 1 Autopilot Connection Diagram. Source [12]. 


• 14 PWM / Servo outputs (8 with failsafe and manual override, 6 auxiliary, high-power 
eompatible) 

• Abundant eonneetivity options for additional peripherals (UART, I2C, CAN) 

• Integrated baekup system for in-flight reeovery and manual override with dedieated 
proeessor and stand-alone power supply (fixed-wing use) 

• Baekup system integrates mixing, providing eonsistent autopilot and manual override 
mixing modes (fixed wing use) 

• Redundant power supply inputs and automatie failover 
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• External safety switch 

• Multicolor LED main visual indicator 

• High-power, multi-tone piezo audio indicator 

• microSD card for high-rate logging over extended periods of time 

4.1.2 Specifications 

The Pixhawk 1 Autopilot has the following specifications as found on [12]: 

Processor 

• 32bit STM32E427 Cortex M4 core with EPU 

• 168 MHz 

• 256 KB RAM 

• 2 MB Plash 

• 32 bit STM32P103 failsafe co-processor 

Sensors 

• ST Micro L3GD20H 16 bit gyroscope 

• ST Micro LSM303D 14 bit accelerometer / magnetometer 

• Invensense MPU 6000 3-axis accelerometer/gyroscope 

• MEAS MS5611 barometer 

Interfaces 

• 5x UART (serial ports), one high-power capable, 2x with HW flow control 

• 2x CAN (one with internal 3.3V transceiver, one on expansion connector) 

• Spektrum DSM / DSM2 / DSM-X® Satellite compatible input 

• Putaba S.BUS® compatible input and output 

• PPM sum signal input 

• RSSI (PWM or voltage) input 

• I2C 

• SPI 

• 3.3 and 6.6V ADC inputs 

• Internal microUSB port and external microUSB port extension 
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4.2 Ground Control Station 

The GCS used for this research was MAVproxy [13]. It is an open-source python based 
GCS which provides flexible communication and command with any autopilot utilizing 
the MAVlink protocol [14]. Even though MAVproxy is written in Python (operating 
system (OS) agnostics language), it was found to be cumbersome to operate the GCS on 
any other platform other than Linux. This is primarily because the source code updates 
quite rapidly to support new features and the core developer (Andrew Tridgell) exclusively 
utilizes MAVproxy in Linux. A significant amount of external libraries are utilized which 
results in a moderate amount of compatibility debugging for other OS’s if desired. 

4.2.1 MAVproxy Features 

The following are summaries which proved to be extremely useful for this research [15]: 

• command-line, console based application. Plugins included in MAVProxy provide a 
basic graphical user interface (GUI). 

• network capable and run over any number of computers. 

• portable; capable of running on any POSIX OS with Python, pyserial, and select() 
function calls, which means Linux, OS X, Windows, and others. 

• light-weight design; runs on small netbooks. 

• tab-completion of commands. 

4.3 Simulation 

The APM environment offers three versions of SITE simulations. The lowest fidelity 
SITE is provided by MAVproxy, which is a simple 6-degree of freedom kinematics model 
with no environment or actuator modeling. This proved to be adequate for initial testing but 
resulted in poorly tuned algorithms when actual flight tests were conducted. The MAVproxy 
simulator was used for basic code debugging but nothing else. 

The second SITE offered in the APM environment is X-plane 10. This is a much higher 
fidelity simulation, which includes actuator models and environmental modeling. X-plane 
10 is primarily used for simulating full-scale aircraft and therefore is difficult to find models, 
which accurately represent the dynamics of small fixed-wing UAS. The open-source 
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community provided model called the “maxi-swift” was similar enough to the airframe in 
this research that it provided adequate SITL modeling which ensured robust flight test. 

The last SITL simulation tested under the APM environment was the RealFlight 7.5 API. 
The RealFlight remote controlled (RC) airplane simulator offers some of the industry’s 
highest fidelity simulations for small aircraft. This product requires an API key to hook into 
MAVproxy. This capability is not yet on the market as of the time of this research, but the 
APM core developers were supportive of this research and ran multiple experiments with 
the RealFlight SITL for early validation. The screenshot in Figure 4.2 was captured while 
conducting testing utilizing the RealFlight SITL. 



Figure 4.2. High Fidelity RealFlight 7.5 Software in the Loop (SITL) 


4.4 Airframe 

The aircraft used for this research was the Flitetest Spear and the Flitetest Explorer [16]. 
The Spear airframe was chosen for its endurance capability of greater than 45 minutes of 
flight time and its large capacity fuselage. The flying-wing architecture keeps the actuation 
requirement to a minimum of two servos by utilizing an elevon configuration. 

Figures 4.3 and 4.4 are example photos from the instructional build website [16]. Figure 4.5 
illustrates the process of building the Spear aircraft. The large blunt nose provides adequate 
space for two 2,200 mAh (12.6volts) lithium polymer batteries wired in parallel. The 
remaining cargo space was used for accommodating the Pixhawk autopilot. 
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Figure 4.3. Spear Airframe. Source [16]. 



Figure 4.4. Spear Cargo Capacity. Source [16]. 


This plane was constructed out of craft foam board. The plans were downloaded from 
flitetest.com [16] and converted to CorelDraw vector files for use in a laser cutter. These 
files were then cut out of four sheets of foam board using the laser cutter. The wing 
halves were joined with standard box tape and hot glue. This provided a cheap and rapid 
construction process which was achievable under four hours of build time. 

4.4.1 Spear Specifications 

• weight without battery: 1.45 lbs (658 g) 

• center of gravity: 3 - 3.5” (76-89 mm) in front of firewall 

• control surface throws: 16°deflection - Expo 30% 

• wingspan: 41 inches (1041 mm) 
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Figure 4.5. Spear Build Process 


• motor: 425 sized, 1200 kv minimum 

• prop: 9 X 4.5 CW (reverse) prop 

• eleetronie speed eontrol (ESC): 30 amp minimum 

• battery: (2) 2200 mAH 12.6 volt minimum 

• servos: (2) 9 gram servos 

The Explorer airframe was ehosen beeause it is a eonventional airframe with highly eoupled 
aerodynamies whieh ean be eonfigured in multiple different failure modes. The sport wing 
provided in the plans was modified to have independently aetuated flaps, and ailerons for 
a eombination of software enabled failure modes. This aireraft was also eonfigured with 
a rudder for testing the lateral aerodynamies eoupling effeets on the adaptive eontroller. 
Eigure 4.6 illustrates the eompleted Explorer airframe used in this researeh. Eigures 4.7 
and 4.8 are examples of the Explorer build proeess and autopilot integration. 


4.4.2 Explorer Specifications 

• weight without battery: 1.08 lbs (493 g) 

• eenter of gravity: 2.25” (57 mm) from leading edge of wing 

• eontrol surfaee throws: 12°defleotion - Expo 30% 

• wingspan: 57 inehes (1447 mm) 

• motor: 425 sized, 1000 kv minimum 

• prop: 9x6 CW (reverse) prop 

• ESC: 30 amp minimum 

• battery: (1) 2200 mAH 12.6 volt minimum 
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servos: (6) 9 gram servos 



Figure 4.6. FliteTest Explorer. Source [16]. 



Figure 4.7. Explorer Cut Foamboard Parts 


The previously discussed COTS autopilot, Linux software tool chain for SITL and GCS, 
and prototype aircraft were utilized to present the results found in Chapter 5. 
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Figure 4.8. Pixhawk Autopilot Installed on Explorer 
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CHAPTER 5: 

Flight Testing and Performance Evaluation 


This chapter utilizes the APM SITL arehiteeture as well as prototype aireraft to tune and 
test the Xi adaptive eontrol algorithm to highlight performanee parameters whieh meet or 
exeeed the underlying requirements established in the problem proposal in Chapter 1. For 
all graphs in Chapter 5, the green baekground identifies when the adaptive eontrol algorithm 
was being utilized. 


5.1 Simulation Results 

Simulation results were eaptured utilizing the APM SITL using X-Plane 10. The aireraft 
model ehosen was an open-souree flying wing model ealled the maxi-swift as seen in 
Figure 5.1, a simulation photo. 



Figure 5.1. Maxi-Swift Flying Wing Model Used in X-PlanelO 


This plane did not afford proper testing of elevon mixing, but the fidelity in the model 
provided utilizing X-Plane 10 was surprisingly aeeurate with respeet to the FT-Spear aireraft 
used for aetual flight test. This drastieally inereased eonfidenee that anything tested in SITL 
would have a high probability of sueeess in aetual flight test. 
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5.1.1 Pitch Attitude Performance 

It can be seen in Figure 5.2 that there exists some steady state error when under PID control. 
This phenomenon is more pronounced when the desired pitch attitude is negative. This 
steady state error is due to the increase in lift caused by the increasing airspeed. The 
PID controller underperforms in this regime. However, the adaptive controller is able to 
compensate for these dynamics fast enough to track the desired attitude with no steady state 
error. The green background in Figure 5.2 identifies when the adaptive control was being 
utilized and the white background was utilizing the PID control. 
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5.1.2 Roll Induced Pitch Disturbance 

When rapidly rolling the maxi-swift aireraft in simulation, there was a notieeable eoupling 
in the piteh axis. The PID eontroller struggles to eorreet this diserepaney beeause the time 
eonstant of the integral error simply eannot be inereased high enough to aehieve satisfaetory 
eompensation. As seen in Figure 5.3, the Xi adaptive eontroller signifieantly redueed the 
pitehing disturbanee due to rapid rolling. 



Figure 5.3. PID vs Xi Adaptive Control Coupled Pitch Response Perfor¬ 
mance 


5.1.3 Pitch due to Landing Gear and Flaps 

Lowering the landing gear and flaps entering the landing phase of flight eauses un- 
eommanded deviation in piteh. The F4U Corsair (see Figure 5.4) in X-Plane 10 was used to 
evaluate the attitude hold retention performanee. This un-modeled aerodynamies ean eause 
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the integrator in a PID controller to saturate or wind up. Figure 5.5 illustrates the pitch 
performance during the actuation of flaps and gear with and with out adaptive control. 



Figure 5.4. F4U Corsair Model in X-PlanelO 


The Xi adaptive controller significantly reduces the attitude excursion due to flaps and 
landing gear. 

5.1.4 Roll Performance 

Adaptive control offered little improvement to roll performance with respect to the nominal 
attitude retention. The roll axis for multiple aircraft was typically seen to have higher cutoff 
frequencies with respect to pitch and therefore required higher values for the D(s) cut off 
frequencies. Figure 5.6 illustrates that the adaptive controller is capable of achieving equal 
or improved performance over the PID controller. 

5.1.5 Yaw to Roll Coupling 

The roll dynamics are coupled to the yaw dynamics as seen in Equation B.l 1. In the case 
of an actuator/servo hardover in the yaw channel, the coupling causes unwanted rolling 
moments. Figure 5.7 compares a left and right yaw servo hard over for both PID and 
the Xi controller. The adaptive controller significantly outperforms the PID controller in 
maintaining the desired roll. It could be argued that the PID controller’s integrator time 
constant could be re-tuned to be comparable. It was evident for each of the simulation 
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Figure 5.5. Pitch Attitude Response Due to Gear and Flaps 


tests that tuning the PID for every seenario would likely have been possible individually to 
aehieve similar performanee to the Xi . However, the Xi was not tuned between eaeh of 
these tests and outperformed PID in most, if not all oases. 

5.1.6 Actuator Miscalibration 

The Xi adaptive eontrol provided fast learning of airframe aotuator misoalibrations. The 
autopilot parameter SERVOX_TRIM was used to offset the ailerons by 15% of their travel 
to evaluate how quiokly the oontroller was oapable of adapting to the new offset. This was 
tested first on the PID oontroller and then on the Xi as seen in Figure 5.8. 

It oan be seen in Figure 5.8 and 5.9 that the new trim value aohieves the referenoe oommand 
in about 0.5 seoonds for the Xi . The Xi is slightly faster than PID, but the Xi exhibits 
some overshoot. It is important to note that this rate of learning is perfeotly adequate 
for learning the misoalibrations in flight, but is insuflioient for learning on take off. In the 
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Figure 5.6. Roll Attitude Performance 

takeoff circumstance, the algorithm saturates the actuator if biases exist between desired and 
achieved. Because the aircraft has no dynamic pressure (it is not flying), the desired cannot 
be achieved. This causes controller saturation, which cannot be unwound fast enough for 
take-off (specifically for tail-dragger configuration). This has to be handled with ad-hoc 
heuristics very delicately. One could choose to speed scale the learned parameters to help 
with learning rate, but it was not chosen for this research because this controller is also 
utilized in tail-sitter configurations where the zero airspeed would cause the controller to 
refuse to learn when in VTOL mode. 


5.2 Flight Test Results 

After successful results were achieved utilizing the SITL environment across multiple 
aircraft, similar test techniques were conducted on the prototype aircraft described in 
Chapter 4. 
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Rudder Servo Hardover/Fallure 
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Figure 5.7. Roll Response to Rudder Servo Hardover/Failure 

5.2.1 General Observations 

The FT Spear aireraft was used in all initial flight tests to prove the algorithm’s initial 
robustness and faeilitate testing of various tuning seenarios. Choosing this aireraft proved 
to be benefieial beeause it provided 30 or more minutes of enduranee. The flying wing 
arehiteeture had a slight forward CG when eonfigured with two 2200 milli amp hour (mAh) 
batteries. This resulted in the aireraft being anemie in piteh response but still adequate for 
testing the Xi algorithm. Various iterations of the algorithm were tested and eompared to 
the results observed in SITL as the development progressed. 

The first observation was that the algorithm, even when poorly tuned, was bounded in 
response and never approaehed instability. The gains found in SITL when applied to aetual 
flight test always resulted in stable flight. As diseussed in Seetion 6.1, the first flight 
tests of gains found in SITL either produeed low-frequeney oseillations or high-frequeney 
oseillations. Low-frequeney oseillations oeeurred when the Xi filter eutoff bandwidth was 
set too low. High-frequeney oseillations oeeurred when the Xi filter eutoff bandwidth was 
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Miscalibrated Aileron Actuator 




Figure 5.9. Xi Fast Adaptation to Unknown Miscalibration 


set too high. Tuning the filter to achieve adequate response somewhere between the low and 
high-frequency oscillations was easily achievable within three to four iterations of setting 
the filter cutolf frequency. Finding the optimal filter setting is arguably difficult to see on 
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the low refresh rate telemetry, but aehieving similar or better performance than PID was 
easily attained with very minimal effort. 

Once the algorithm was well understood and properly implemented in code, flight tests were 
then conducted on the FT Explorer aircraft. This aircraft provided more test opportunities 
because multiple aircraft configurations could be tested with various combinations of flaps 
and rudder. 

5.2.2 Pitch Performance 

As demonstrated in Figure 5.10, the simulation results almost identically matched the flight 
test results for pitch response. As the aircraft pitches down, the airspeed builds and causes a 
building steady state error that the integrator time constant struggles to out pace. This was 
observed both in simulation and in flight test. The Xi adaptive controller’s fast adaptation 
was sufficient to counteract this unmodeled pitching moment. These results drastically 
increased the confidence both in the SITE fidelity as well as the Xi adaptive controller. 

5.2.3 Pitch due to Landing Gear and Flaps 

The Xi adaptive controller performed well when compensating for the change in aircraft 
configuration when lowering the flaps. Figure 5.11 illustrates that the tracking error from 
lowering and raising the flaps under adaptive control is almost imperceptible. However, 
PID resulted in a sharp peak both when lowering and raising the flaps. 

5.2.4 Roll Induced Pitch Disturbance 

The aircraft max roll angle was set to +45° to achieve rapid roll while trying to maintain a 
zero pitch attitude. This rapid roll maneuver caused significant excursions in pitch when 
under PID. However, the Xi controller maintains the pitch attitude within ±3°. It can be 
noticed in Figure 5.12, that the adaptive controller oscillates around zero pitch attitude, 
while the PID has more random excursions. The slight sinusoidal behavior or the Xi 
controller is likely due to the cutoff frequency of the filter being slightly too low. This 
maneuver may possibly be utilized to fine tune the Xi’s filter cutoff frequency for pitch. 
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5.2.5 Roll Performance 

To test the nominal flight path stability and roll attitude performanee, a reetangular flight plan 
was established. Figure 5.13 illustrates the Xi adaptive controller while maintaining the left 
hand pattern. The noise in the roll channel was observed in actual flight as turbulence. The 
filter cutoff frequency could be lowered to reduce performance if desired, but these results 
were not oscillatory which suggested that the filter cutoff frequency was not too high. These 
results compare to the noise found on the PID channel so the threat of damaging servo by 
over driving them was not expected. 

Asymmetric flaps were configured by only allowing one servo to lower the left flap while 
leaving the right flap stationary. This failure was designed to test the unmodeled roll 
response from a very asymmetric mechanical failure. Unfortunately, the roll and pitch 
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Pitch Performance Deploying Flaps with and without Adaptive Control 



excursions caused by this failure were marginal. Both PID and the Xi handled this failure 
with adequate reliability. 

The overall result of the performance testing of the Xi algorithm was quite successful. The 
results from SITL provided confidence in the expected results from actual flight test and 
proved to be robust representations of the results achieved in flight test. Flight test were 
significantly more difficult to conduct, but the mirrored performance between the SITL 
and flight tests proved that more time can be spent utilizing SITL to drastically reduce 
development cost and schedule. 
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Figure 5.12. PID vs Xi Adaptive Control Coupled Pitch Response Perfor¬ 
mance 
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Figure 5.13. Xi Adaptive Control Roll Performance 
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CHAPTER 6: 
Recommendation 


This chapter provides various recommendations pertinent to the integration of the Xi adap¬ 
tive control algorithm. These recommendations may be more qualitative than quantitative, 
but provide experience based guidance for the engineer to better understand the nuances of 
the Xi adaptive control algorithm. Prior to this research, there was very little experientially 
based guidance on implementation of this algorithm which proved to be the largest barrier 
for implementing this modern controller. 

6.1 £i Adaptive Control Algorithm Tuning 

The primary goal of this research was to reduce the complexity of tuning expertise required 
to get an unknown airframe airborne successfully. The amount of time required to get 
the adequate airframe performance using the Xi algorithm is significantly reduced. Three 
primary features need tuning: the adaptation gain, the controller filter, and the companion 
model. 

6.1.1 lYining the Adaptation gain 

Tuning the adaptation gain is fairly intuitive. The adaptation gain is a function of the loop 
frequency at which the filter is run so it may only need to be tuned for a given autopilot 
at a given loop rate and not for every specific aircraft. Anecdotally, the adaptation gain 
of 10,000 was used on the Pixhawkl autopilot running the scheduled loop rate at lOOHz 
across multiple aircraft without needing to modify. The primary feedback to the user for 
tuning this gain resides in the desired movement of the estimated states. The adaptation 
gain was set low (1-10) while watching the parameters (6, a>, cr) adapt. The adaptation gain 
is slowly increased until the desired rate of adaptation as seen in the real time monitoring of 
the parameters is adequate. Another approach for tuning the adaptation rate is to increase 
the gain until parameters start to oscillate between the bounds and then reduce it. This 
bang-bang response in the parameter adaptation will show up in the performance of the 
controller as increased peaks away from zero error between desired and achieved. In the 
case of this research, this noise can be hard to identify in the rate control itself and therefore 
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is why monitoring the attitude error helps find the appropriate adaptation gain whieh is 
extremely high but still not injeeting attitude error spikes. 

6.1.2 lYining the Controller Filter 

The adaptation feedbaek gain (k), as diseussed in Seetion 3.3 was by far the most influential 
gain to tune. The simplieity of tuning the Xi algorithm resides in the faet that the majority of 
lay users eould adequately tune an airframe with this stand-alone gain. Adequate values for 
this gain ranged from 0.3 for responsive aireraft and 2.0 for very sluggish aireraft. This value 
was set for both the roll and piteh axis independently. As seen in Equation 3.8, this value 
is establishing the eutoff frequeney of the eontrol filter that separates the bandwidth limited 
eontrol ehannel output and the high bandwidth adaptation. The default value assigned in 
the souree eode was 0.45, whieh proved to be a good starting point. If the default value 
for k is not eorreet, then the eontrol ehannel will produee either low-frequeney oseillations 
or high-frequeney oseillations. The low-frequeney oseillations are produeed beeause the 
bandwidth of the eontrol is too low and there should be a pereeptible lag between the 
desired state and the aehieved state. High-frequeney oseillations oeeur when the eontrol 
filter bandwidth is set higher than the plant’s bandwidth, and the aireraft is ineapable of 
aehieving the desired rates. Unlike PID eontrol, the Xi eontrol never exhibited unstable 
performanee with ineorreet gains. The eontroller simply oseillates with extremely poor 
performanee. This was a remarkable feature beeause poorly tuned PID gains ean eause an 
aireraft to depart eontrolled flight rapidly. Whereas, the Xi eontroller maintained bounded 
flight performanee as the theory suggests. 

6.1.3 Timing the Companion Model Cutoff Frequency 

System identifieation was eondueted as seen in Appendix C in order to aseertain the 
bandwidth of various airframes and their aetuators. Figure 6.1 is an example of seeond- 
order models for two aireraft’s roll dynamies eompared to seeond-order models of RC 
aetuators (servos). As to be expeeted, the bandwidth for the aetuators is slightly higher than 
the airframe dynamies. 

These rough approximations were used to then plaee the eutoff frequeneies of the eompanion 
model with the expeetation that the eompanion model eannot aehieve higher bandwidths 
than that of the airframe. Conservatively, the eompanion model eutoff frequeney was 
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typically set 2-5 rad/s lower than the expeeted max performanee of the airframe. After the 
other algorithm gains are tuned, and satisfaetory performanee is aehieved, the eompanion 
model eutoff frequeney ean then be inereased to aehieve higher performanee. 


6.2 Improved Recursive Architecture 

The speed of adaptation and aeeuraey of the diseretized Xi algorithm is drastieally improved 
with inereased loop frequeney. The algorithm was written as one reeursive arehiteeture that 
updates at the Pixhawk’s seheduled loop rate. As previously diseussed, the seheduled loop 
rate ideally should run at lower frequeneies to prevent exeessive log file size and added strain 
on the CPU. However, the Xi arehiteeture only requires the adaptation loop be run at faster 
rates. With this speeifie performanee enhaneement in mind, the APM arehiteeture eould 
be modified to aeeept an independent loop speeifioally for the Xi adaptation update. The 
sensors measurements and extended Kalman filter (EKF) update ean run signifieantly lower 
with only inereased performanee of the algorithm. The higher adaptation loop would enable 
higher adaptation gains (T) and eonsequently produee faster adaptation of the system. 
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6.3 Integrator Windup Issue 

The Xi controller architecture exhibits a similar response to integrator wind up in PID 
control when the aircraft is not flying. If the controller was enabled before takeoff, the control 
surfaces would move to counteract any slight deviation from input to output attempting to 
zero the error. Because the aircraft was not flying, the controller would continue to increase 
the control output until the actuators reached saturation. The baseline code was configured 
to ensure that the parameter estimates would not continue to integrate after saturation was 
reached. This technique is standard practice when writing control laws that utilize actuators 
that have saturation limitations. However, the anti-windup feature that this offers only 
occurs when the actuators are saturated. This is inadequate for the takeoff scenario. The 
flight surfaces quickly saturate while waiting for takeoff and cause a crash immediately 
upon takeoff. It is standard practice to also disable the integrator action of controllers if it 
is known that the aircraft is not flying using any combination of airspeed, throttle position, 
inertial measurement unit (IMU) estimated velocity, or global positioning system (GPS) 
ground speed. In the case of the Xi algorithm, the integrator utilized in D(s) is essential 
to the entire architecture and presented challenges in how to use ad-hoc methods to prevent 
the integrator windup issue even though the “non-flying” state was calculable onboard the 
autopilot. The initial experiments were to see if the adaptation rate was fast enough to 
un-learn the aircraft’s saturated state, but there were no combinations of filter gains which 
resulted in learning rates adequate to ensure safe takeoffs. Further research in this area is 
required to completely replace the PID controller with the Xi . All takeoffs were conducted 
either in manual control or with the PID controller enabled until safely airborne. 

In summary, no one manual has been created for this type of controller’s implementation and 
integration. The recommendations provided may only apply to this specific implementation 
of the Xi adaptive controller, but it offers guidance where none previously was articulated 
in contemporary literature. 
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CHAPTER 7: 
Conclusion 


The primary objective of this thesis was to evaluate a modern control technique to determine 
if the advanced controller could reduce the Navy’s growing GNC demand on engineering 
resources. 

The Xi adaptive control algorithm suggested and developed in this thesis proved to be a 
successful alternative to the conventional PID control strategy which drastically reduced the 
cost and time requirements to achieve robust flight. In addition to meeting the objective of 
lowering the engineering demand of the aircraft GNC algorithms, the Xi adaptive control 
could also provide battle damage tolerance and achieve more accurate control with the 
utilization of fast and robust adaptation. 

Conversion of the continuous time Xi algorithm proved to be difficult with limited precedent 
literature but proved to be quite achievable with very little limitation posed by contemporary 
embedded processors. The APM flight stack open source project proved to be a very versatile 
software base that enable rapid prototyping and testing of the Xi algorithm. 

The use of SITL was a critical tool in the development of this algorithm into source 
code. Follow on research or implementation of aircraft GNC in the simulation environment 
drastically reduces the project risk at little to no expense using contemporary high fidelity 
aircraft simulation tools. 

Flight test of the Xi algorithm highlighted many capabilities gained through the use of 
fast and robust adaptation. The primary goal of reducing the GNC engineering demand 
required to achieve successful robust flight of UAS’s was evident even on the first test 
flights. Achieving adequate flight performance was now attainable within a matter of 
minutes instead of multiple tests spread across multiple flights potentially spread across 
multiple days. With the use of SITL integration, a robust model of the airframe, and more 
dedicated research, successful first flights with no tuning required are now within the realm 
of possibility. 
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APPENDIX A: 
Transfer Functions 


A.l Transfer Functions 

This research utilizes the transfer function (TF) representation of aircraft flight dynamics 
that is typical of linear time invariant (LTI) systems. A transfer function is a very useful 
approach to describe the relationship between inputs and outputs of LTI systems. Both 
analytically and numerically, the TF approach has significant benefits in continuous and 
discrete time domains as its construct is based on well-developed properties and primitives 
of polynomials. These polynomial representations in the s or z domain map to aerodynamic 
and inertial coefficients through equations of motion. What is unknown or partially known 
a priori, are the numerical values of coefficients for those polynomials. Therefore the tools 
from the areas of online estimation such as regression are utilized to solve for them. 

Transfer functions take the form 


H{s) = 


TO 


(A.l) 


where 


Y(s) is the Laplace transform of the output 
X(s) is the Laplace transform of the input 

Standard physics models of first and second order form are well understood and seen in 
many model derivations. The first-order model takes the form 

H(s) = (A.2) 

T5 -I- 1 S + CO„ 

where 

kcic is the DC gain 
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T is the system time constant (time in seconds to reach 63% of steady state) 
a)„ is the natural frequency of the first order system 

Similarly the standard form for a second-order system takes the form 


H{s) = 


COn 


+ I^COqS + Cl>\ 


(A.3) 


where 


a>o is the system natural frequency in radians per second 
4" is the system damping ratio 

The modeling of a system can also be converted to a system of first-order differential 
equations also known as state-space modeling. In this case, the first-order system model 
can be represented as 

x{t) = Ax{t) + Bu(t) (A.4) 

where x is the time derivative of the state, A is the state transition matrix with all its 
eigenvalues chosen negative (Hurwitz), B is the input matrix, and u is the input vector. 


60 



APPENDIX B: 

Fixed Wing Aircraft Dynamics Model 


B.l Fixed Wing Aircraft Dynamics Model 

The following is the nomenelature used to deseribe the fixed-wing kinematie equations. 
Euler angles for piteh (9), roll (0), and yaw (0) have the units of radians. Figure B.l 
illustrates the north-east-down (NED) referenee frame definitions used for body rotational 
rates about the x axis (p), y axis (q), and the z axis (r), as well as the body veloeities in 
the X axis (u), y axis (v), and the z axis (w). 



Figure B.l. Reference Frame - Body Rates and Velocities. Adapted from 
[16]. 

Newton’s seeond law as it pertains to rotational motion ean be stated as 


where t is the torques applied to the body, J is the moment of inertia, and ^ is the angular 
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acceleration of the body with respeet to the inertial frame. 

Equation B.l can be rewritten in the body referenee frame as follows: 

r‘= X (B.2) 

The expression is the angular aeeeleration in the body frame as viewed in the body 
frame 


4,i = 


‘p^ 




The equation ean then be rewritten with respeet to in Equation B.3 


(B.3) 



= 7-1 


-to 


bfi 


X 


(H-) 


+ T" 


where 7 is the inertia matrix as follows: 


(B.4) 


Jx 

-Jxy 

Jxz 

-Jxy 

Jy 

-Jyz 

~Jxz 

-Jyz 

Jz 


(B.5) 


The moments of inertia, or the diagonal terms, are always non-zero for any rigid body. The 
products of inertia, or the off-diagonal terms, are terms which describe the inertial eoupling 
between axis. Eor a traditional fixed wing aireraft, the natural symmetry will simplify the 
inertia matrix in the off-diagonal terms as follows: 


7 = 


/ 




Jx 

0 


0 -J,z 

Jy 0 
0 Jz ^ 


(B.6) 
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The inverse of J can be found to be 


o 

u p 

T 0 (B.7) 

0 — 

where 

r = 7,7, - Jl (B.8) 

Aircraft nomenclature for torques are defined t = (/, m, n)^, and therefore the combined 
equations derived from first principles take the form 





T 0 W J,,pq + (Jy - J,)qr I 

0 0 Jxz(r^ - P^) + Uz- Jx)pr + m 

^ 0 ^11 (Jx- Jy)pq - Jxz^r U 


Tipq - Y2qr + r3/ + r4n 
Yspr - Y(,(p^ - r^) + j-m 
Y-ipq - Y\qr + r4/ + Ygn 


where 


(B.9) 
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r 


(B.IO) 


The aerodynamic torques (excluding propulsive torques) can be found to be 



= \py;s 

b [C/g + Cl^/3 + C/p 2VaP 2Va ^ + Clsr ^ r] 


m 

C [^Cmg "I" Cnia^ 

(B.ll) 



\b \Cno + CnpP + C„p + ^'J<Sa br\ j 



where C/^, C/^, C/^, C/^, C/^^, C/^^, C„o, C„^, C„^, are the lateral aerodynamic 

coefficients 

and Cmg, Cm„, are the longitudinal aerodynamic coefficients. 

Substituting the aerodynamic torques found in Equation B.ll into Equation B.9 results 
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in [17] 


P — r2^r + 2^^! ^PQ + + ^Pp 2y '*' ^Pr 2y ^PSa "*" ^PSr 

9 9 1 9 1 

4 ~ ~ ^6iP ~ ^ ) ~^ '^P^a ^^0 "*” Crua^ "*" ^tT *” 

r = Tvi?^ - Tiqr + + Cr^/3 + Crp^ + + Cr^,5r 

(B.12) 

Simplifying Equation B.12 assuming no inertial or aerodynamic coupling results in 

P = ^PSa ■*■ ^Pp 2^ 

4 = (B-13) 

^ = 2 ^fSr ^r + C^r ■*■ ^ro 

The equations in B. 13 are then slightly modified to fit the first-order ordinary differential 
equation (ODE) model as described in Equation A.4 

p = App -I- bp [d)p6a + Opp + d-p') 

q = Aqq + bqicjqbe + bqq + o-q) (B.14) 

f = A^r -I- bp (^d)r6r + Of-r + d-p^ 

where to is the input gain coefficient, 9 is the constant state coefficient, and cr is the 
disturbance estimate. 
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APPENDIX C: 
System Identification 


C.l System Identification 

System identification was performed to scope the author’s expectation for the bandwidth of 
small fixed wing unmanned airframes. This was conducted in order to properly establish 
cutoff frequencies for the companion model. These results are also mirrored in Figure 6.1 
and discussed in Chapter 6. 


C.1.1 Data Collection 

The Pixhawk autopilot was used to capture roll and pitch rates (p, q) for the test vehicle as 
well as the pilot’s command inputs. These outputs and inputs were the essential building 
blocks for creating pitch rate and roll rate models for the test vehicle. The autopilot is 
capable of logging data at 50-400 Hz that is represented as a discrete time domain signal. 
This data should ultimately be manipulated into the s-domain. The mathematics for this 
procedure are well defined, and numerous tools can be used to simplify this process [18]. 

It is crucial to ensure there is sufficient frequency content in the data recorded. Exciting 
multiple frequencies in the time domain ensures the regression techniques have an adequate 
persistence of excitation to resolve polynomial coefficients with higher certainty. 

To ensure sufficient frequency content was obtained from the aircraft, a linear chirp was 
chosen and implemented into the Pixhawk source code as follows: 


x{t) = sin 


(pQ-\-ln\ fot + 


(C.l) 


where 


00 is the initial phase of the chirp at t=0 (nominally zero degrees) 
/o is the initial frequency at t=0 
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k is the chirp rate 
t is time in seconds 

An example of this method can be seen in Figure C.l. 



Figure C.l. Reverse Linear Chirp Example 


Specifically for this research, the reverse formulation of the linear chirp was used because 
it ensured the aircraft wouldn’t exceed max roll or pitch limits within the first few cycles 
of the output. The reverse chirp was applied in open loop while in “training mode” which 
would constrain the angle of bank or pitch by switching into attitude hold mode upon an 
attitude exceedance. 

C.1.2 z-Domain to s-Domain 

The logged input and output data in discrete form requires shaping to convert cleanly into 
an s-domain representation. The first step is ensuring that the data is of constant sampling 
rate. In other words, the time between samples is uniform from sample to sample. The data 
provided from the Pixhawk autopilot does not have a uniform sampling rate. The sample 
rate is a user-defined rate (50-400Hz) but has a slight amount of jitter (+0.1%). piecewise 
cubic hermite interpolating polynomial (PCHIP) interpolation was used to interpolate the 
data into a uniform sampling rate. 

After the data is shaped correctly with a uniform sample rate, the discrete data is transformed 
into a continuous approximation using a zero order hold (ZOH) technique. Taking the 
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Laplace transform of the continuous input/output data will convert it into the s-domain, and 
finally, a non-linear least squares minimization algorithm can be run to find the polynomial 
coefficients which best fit the data. 

The order of the regression (number of polynomials to estimate) is at the discretion of 
the engineer and their intuition of system’s physical representation. Higher order models 
will better fit the data, but in most cases, they tend to overfit the data; therefore some 
tradeoffs should be considered to simplify the model. Most aircraft models reasonably limit 
the system to an LTI and second-order. These fundamental aerodynamics models divide 
the modeling into longitudinal and lateral dynamics. Each axis of the aircraft is assigned 
two-second order responses. Pitch, for example, has a second order response in the pitch 
damping mode (also known as the short period) and also has a second order response in the 
transition of kinetic energy to potential energy (also known as the long-period or phugoid). 
Both first order and second order models were estimated for comparison sake. 

Results were collected from two flight test events. The first flight test was conducted under 
manual control without utilizing the chirp. The pilot attempted to increase frequency of 
the input signal manually. The second set of data collected was via the reverse linear chirp 
method previously described. 

The manually piloted acquired data was expected to have insufficient frequency content 
in the signal. However, the manual flight test provided moderate results for modeling the 
aircraft as seen in Figure C.2 . 

The results in Figure C.2 demonstrates the utility of this technique even if data can only 
be acquired from manual pilot inputs. It can be seen that the second order model starts 
to misrepresent the data at higher frequencies. Figure C.3 illustrates the data recorded 
from the reverse chirp experiment run at 50Hz. The S pulse width modulation (PWM) input 
channel data clearly does not represent the software calculated chirp. The highest frequency 
designed to be outputted during this experiment was lOHz. 10 Hz is significantly less than 
the frequency required to meet Nyquist sampling criterion (25Hz) for this data sampling 
rate. However, there are clear patterns of aliasing in the input signal as recorded by the data 
flash logger. 

The chirp response was physically observed on pre-flight, in actual flight, and in the data- 
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Roll Analysis 




Figure C.2. Roll Model Regression with Manual Inputs 




Figure C.3. Roll Model Regression with Reverse Linear Chirp 
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flash logged body rate of the aircraft. However, the aliased input channel was arbitrarily 
biasing the regression result. The significant aliasing in the logged input was not due to 
insufficient sampling rate as previously discussed. After inspection, the peculiar aliasing 
issue was hardware specific to the Pixhawk 1 autopilot in how the main CPU sends servo 
commands to the auxiliary I/O CPU. The most recent version of firmware, at the time of 
this test, improperly logs the PWM through an aliased prone signal path. The main and I/O 
CPU both run at 50Hz with some appreciable clock drift. This generates a noticeable beat 
frequency and delay when the actual values in registry are saved for PWM values are sent 
back round trip to the main CPU. The implication of logging the PWM values at the very 
end of the digital transmission line seems valuable in principle because the values being 
logged are the undeniable values being sent to the actuators. However, the cost of logging 
these values in this manner on the Pixhawk architecture incurs significant aliasing at almost 
all frequencies. Logging the commanded PWM values before being sent to I/O CPU solved 
the aliasing discrepancy and produced very frequency rich models. 

The manually piloted acquired data provided viable data source for the models even though 
it is a very simplistic approach. There were two separate manual tests run on the same 
aircraft on the same flight, and the following are the results using this regression technique 
to model a second-order system: 


H(s) = 


10.39 

52 + 31.265 + 504.9 


(C.2) 


and 


H(s) = 


10.61 

52 + 29.775 + 498.7 


Converting to standard form as described in Equation A.3 yields 


H(s) = 


0.0206 * 22.472 

52 + 2 * 0.69 * 22.475 + 22.472 


(C.3) 


(C.4) 


and 
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H(s) = 


0.0213 * 22.332 
s^ + 2* 0.67 * 22.335 + 22.332 


(C.5) 


It is important to note that this system identification technique run on separate sets of data 
has produced two models with very similar values for and 

This produces the average values of 

a)„ = llArad/s 
k = 0.0209 
4" = 0.681 


With the aliasing removed from the chirped input command signals as previously described, 
the model is drastically improved and produces the following results: 


H(s) = 


4.409 

52 + 27.115+ 430.6 


(C.6) 


and 


H(s) = 


3.295 


f2 + 18.825 + 296.5 


(C.7) 


Converting to standard form as described in Equation A. 3 results in 


H(s) 


0.0102 * 20.752 


+ 2 * 0.65 * 20.755 + 20.752 


(C.8) 


and 


H(s) = 


0.0111 * 17.2E 


?2+ 2*0.54* 17.215 + 17.212 


(C.9) 


The comparison of the model to recorded results is seen in Figure C.4. 


72 





Figure C.4. Non-Aliased Reverse Chirp Model Example 


This produces the average values of 

iOfi = 18.98raJ/j' 
k = 0.010 
^ = 0.598 


In the author’s experience, these values are reasonable values for this size and weight 
of airframe based on multiple system identification experiments in SITL across multiple 
airframes as well as actual test flights of three different airframes. This regression technique 
has shown potential to create realistie models from actual flight test data. The data must be 
properly shaped. The reverse ehirp method has the potential to increase the fidelity of the 
high-frequency response of the aireraft if the aliasing issue can be resolved on the eommand 
input ehannel. 
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APPENDIX D; 
Lyapunov Stability Definition 


Aerospace designs tend to use linear controllers as reference models for their simplicity and 
well-understood robustness characteristics. This is despite the fact that the applications of 
these linear controllers are applied to a non-linear dynamical system such as attitude control 
with varying dynamic pressure. Adaptive Controllers are non-linear and may offer per¬ 
formance benefits to non-linear systems as seen in aforementioned aerospace applications. 
However, non-linear controller’s stability properties need to be evaluated. The Lyapunov 
stability criteria offer methods of evaluating these controller’s boundedness and robustness 
behavior. 

Aleksandr Lyapunov was a Russian mathematician who’s work was published in 1892 [19] 
concerning the behavior of non-linear systems close to equilibrium without having to 
rigorously find the unique solutions to difficult differential equations used to model the 
system. His work was largely overlooked until the Cold War when aerospace solutions 
required a more rigorous approach to analyzing non-linear control robustness. Modern non¬ 
linear control engineers extensively utilize Lyapunov’s techniques to design and evaluate 
non-linear controllers. 


D.l Lyapunov Stability Theory 

According to Lyapunov, the stability properties of a system can be classified as stable, 
asymptotically stable, and exponentially stable. 

Given the following autonomous nonlinear system: 


x{t) = f(x(t)), x(0) = XQ 


(D.l) 


where / has equilibrium at Xe : 


f(Xe) = 0 


(D.2) 
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then the equilibrium is said to be: 

1. Lyapunov Stable 

for every e > 0 there exists a d > 0 such that, if ||x(0) - Xe\\ < 6, then for every 
t > Owe have||v(0 - Xe\\ < e 

2. Asymptotically Stable 

if the system is Lyapunov stable and there exists a d > 0 such that if ||x(0) - Xe\\ < 6, 
then lim||x(0 - XgW = 0 

t—^oo 

3. Exponentially Stable 

if the system is asymptotically stable and there exists a>0, /3>0,6>0 such that if 
||x(0) - XeW < d, then \\x(t) - Xe\\ < crllxCO) - Xe\\ for all i > 0 

Being Lyapunov stable infers that if a system is near equilibrium, then it will indefinitely 
remain near equilibrium. If the system if found to be asymptotically stable then it eventu¬ 
ally will achieve equilibrium as i ^ oo and being exponentially stable implies it reaches 
equilibrium with a rate of convergence /S. 

Lyapunov’s Second Method 

Lyapunov’s second method is also known as Lyapunov stability criteria. This method offers 
a less tenuous method for evaluating mathematically non-ideal systems. Lyapunov analysis 
of the linearized system around equilibrium can be cumbersome when the eigenvalues are 
purely imaginary. In this case, the solutions can rapidly depart to infinity or approach zero 
with little perturbation to the eigenvalues. Lyapunov’s second method offers an alternative 
approach for classifying a system’s stability using a concept that is similar to how total 
energy is defined in classical mechanics. 

Conceptually, Lyapunov’s second method can be compared to the evaluation of mechanical 
system similar to the modeling of energy in a vibrating spring mass system. The energy of 
the unforced spring mass system will dissipate energy due to friction and or damping. This 
trend of energy leaving the system towards some “attractor” is evidence of the system’s 
stability characteristics and identifies that there will be some stable end state. Likewise, 
Lyapunov’s second method specifies this with the use of a Lyapunov candidate function 
V (x) which implicitly characterizes the total energy of the system. It is important to note 
that Lyapunov realized that the candidate function could be any positive definite and radially 
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unbounded function. It is then said to be Lyapunov stable if any candidate function is found 
and meets the stability criteria. This means that it is only incumbent upon the engineer to 
find one candidate function to meet the criteria. This can be an iterative process of trying 
various energy like equations. A common approach is to model the Lyapunov candidate 
equation as kinetic energy Lyapunov realized that characterizing the energy of a 

nonlinear system could be almost impossible for some cases, but this approach could prove 
stability without the rigorous knowledge of the true system’s energy. 

Lyapunov’s second method defines a system as Lyapunov Stable for a system i = f(x) 
having an equilibrium point at .r = 0 where the Lyapunov candidate function y(v:) : R" ^ R 
such that: 

• y (v) = 0 if and only if .r = 0 

• y(.r) > 0 if and only if x 0 

• V(x) = ^(V(x) = 2 ^ 0, for all values of .r 9^= 0 

;=1 ' 

if V(x) < 0 for jc 0 then system is asymptotically stable. 

Additionally, it is required to demonstrate the condition of radial unboundedness to ensure 
the system is globally stable [20]. 


77 



THIS PAGE INTENTIONALLY LEET BLANK 


78 



APPENDIX E: 
Projection Operator 


E.l Derivation of Projection Operator Proj{6, y) 

In this section, the projection operator will be defined. Adaptive controllers often use the 
projection operator in their adaptive laws to ensure uniform boundedness of the system error. 
This can aid in faster adaptation and ensure controller closed-loop stability. The projection 
operator attempts to mathematically achieve two objectives; ensure the Lyapunov function 
time derivative remains negative semi-definite, and to keep the estimated parameters uni¬ 
formly bounded. Using the projection operator ensures that 6 is locally Lipchitz continuous 
even though the input y is piecewise continuous. The projection operator as utilized in this 
research is defined as 


Projie, y) = 


y 

' y 

V/ 

^ IIV/II 


(iiv/ir 


> 0 , 

if f (6) < 0 and Vf^y > 0, 
if f (6) < 0 and Vf^y < 0. 


(E.l) 


where 6 > 0 


f(0) = 


02 


-e 


2 

max 


60 


2 

max 


vr = 


20 


60 


2 

max 


(E.2) 

(E.3) 


The projection function chosen for this research is parabolic and has inflection points at user 
defined maximum/minimum bands. Equation E.2 is plotted in Eigure E.l with example 
maximum/minimum of 0.65 and various values for e. The engineer must set the value 
for 6 to achieve the desired slope of the projection at the maximum values. This slope 
should be steep enough to capture the highest expected error given one recursion through 
the algorithm. 
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e 


Figure E.l. Projection Operator - f{9) 


There may exists systems whieh need to have projeetion bounds whieh are not symmetrieal 
about zero. This was found to be the ease for this researeh and therefore the following 
projeetion funetion was used to assist the engineer tuning the algorithm: 


where e > 0 


f{e) = - 

.T 4(6, 


4(6niin - 


0) 


V/' 


^(^max ^min) 

+ ^max ~ 26 ) 
e(0max - ^min)^ 


(E.4) 

(E.5) 


Equation E.4 is plotted in Eigure E.2 with example maximum of 0.65, minimum of 0.25, 
and various values for e. 
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Projection Operator 



Figure E.2. Projection Operator with Offset Limits 


E.2 C++ Implementation 

float proj e c t ion_oper a t or ( f loa t theta , float y, float epsilon, 
float theta_max , float theta_min) 

{ 


// Calculate convex function 

float f_theta = (-4*( theta_min-theta )*( theta_max-theta )) 
/( epsilon *(theta_max-theta_min 

float f_theta_dot = (4*( theta_min + theta_max-(2* theta ))) 
/( epsilon *(theta_max-theta_min 

float proj ec ti on_ou t = y; 

if (f_theta <=0 && f_theta_dot*y < 0) 

{ 
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projection_out = y *( f_t he t a + 1); // y-(y*-f_t he t a ) ; 


} 

return projection_out ; 
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APPENDIX F; 
Matlab Code 


F.l Euler vs Trapezoidal Method 

clear , format long , clc , close all 
dt = 1; 

to = [0:.01;4]; 

yO = exp (to ) ; 

%E uler integration 

tl = [0:4]; 
yl = exp (0) ; 
for i = 1:4 

yl ( i+l)=yl (i) + dt*exp(i-l); 

end 

%Trapezoidal integration 
t2 = [0:4]; 
y2 = exp (0) ; 

for i = 1:4 

y2( i+I)=y2( i) + dt*exp(i-l); 
y2(i+l)=y2(i) + (dt/2)*( exp (i -1 )+y2(i+1)); 

end 

plot(t0 ,y0, ’k—■ ,tl ,yl , ’b’ ,t2 ,y2, ’r’) 

legend (’ y=e^ tEuler MethodTrapezoidal Method’) 


F.2 SISO Lyapunov Solution Proof 

clear , clc , format compact, close all 
Wo Find Pb 

wn = [5:0.1:10]; %r ad / s 
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for i=l:length( wn) 
b = [ wn (i)]; 
a = [1, wn (i)]; 

[A,B,C,D] = tf2ss(b,a) 

Pb(i) = lyap(A,l) ; 

end 

Pb_test = l./(2if:wn); 

plot (wn,Pb_test , ’o' ,wn,Pb) 

%in this simple 1x1 matrix case Pb is easy to calc by hand in code 
% Pb ends up being: Pb = l/2H^wn; 


F.3 Reverse Linear Chirp 

clear , format compact , clc , close all 

fs = 1/200; 

t = [0;fs :7]; 

f0= O.OliWz 

fl= 10; %nz 

k = (fl-fO)/(t (end)) ; 

phi_0 = 0.0; 

sample_delay = 0.02; 

out = sin(phi_0+2*pi*(f0*(7-(t + sample_delay))+(k/2) .:i:(7-(t + sample_delay) 

).^ 2 )); 

%out = 1500 + out*10; 
out = out*10; 
plot (t,out) 


F.4 Projection Operator Example Plots 

clear , clc , format compact, close all 
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Plot Projection 

epsilon = [0.25,0.28,0.3,0.38]; 
theta_max = 0.65; 
theta_min = 0.25; 

center = (theta_max + theta_min )/2 ; 

span = (theta_max-theta_min ) 1.4; 

theta = [center-(span/2);0.01;center+(span/2)]; 

for j = 1: length ( epsilon ) 
for i = 1: length (theta) 

%f_the t a(i,J) = -(theta(i).^2-theta_max^2)/(epsilon(J)H=theta_max 
^2) ; 

%f_theta_dot (i , J ) = -(2^^ theta (i)) /( epsilon (j ) theta_max ^2) ; 
f_theta(i,j) = -(theta_min-theta (i)) theta_max-theta (i))./( 
epsilon (j )) ; 

f_theta_dot(i,j) = (theta_min+ theta_max -(2*theta (i))) / ( epsilon (j 

)); 

end 

end 

figure 

for j = 1: length ( epsilon ) 
hold on 

plot(theta ,f_theta(:,j)) 
end 

line ([ min (theta) ,max( theta)] ,[0,0] , ’color’ ,’blue’ ,’linestyle ’ , ’ — ’ ) 
line ([theta_min ,theta_min] ,[ min (f_ theta (: ,1)) ,max( f_ theta (: ,1))] , ’color’ 
,’red’,’linestyle’, ’ — ’ ) 

line ([ theta_max ,theta_max] ,[ min (f_ theta (: ,1)) ,max( f_ theta (: ,1))] , ’color’ 
,’red’,’linestyle’, ’ — ’ ) 

legend(’\epsilon=0.25’ ,’\epsilon=0.28’ ,’\epsilon=0.30’ ,’\epsilon=0.38’) 

ylabel(’f(\ theta)’) 

xlabel(’\ theta’) 

t i 11 e ( ’ Pro] ec ti on Operator’) 

hold off 


F.5 System Identihcation 
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clear , clc , format compact , close all 


% Import data file for analysis 

[filename, pathname] = u i g e t f i 1 e (’ .m’ Choose first MATLAB file'); 
run ( filename ) 

Wo 

% Clean tigger PWMs 

for r = l: length (RCIN. data (: ,9) ) 

if RCIN. data (r ,8) < 1700 %channel 6 
RCIN. data(r ,8) = 1000; 

end 

if RCIN. data (r ,9) < 1700 ^channel 7 
RCIN. data(r ,9) = 1000; 

end 

end 

%Trigger Channel Data 
roll_trigger_data = RCIN . data (; ,9) ; 
pitch_trigger_data = RCIN. data (: , 8) ; 
trigger_data_time = RCIN . data (: , 2) ; 

^Command Data 

rolI_command_data = ((RCOU. data (: ,3) -1500)+(RCOU. data (: ,4) -1500))+1526; 
pitch_command_data = ((RCOU. data (: ,3) -1500)-(RCOU. data (: ,4) -1500))+1460; 
command_data_time = RCOU. data (: ,2) ; 

WMU Data 

p = IMU2. data (; ,3) ; 
q = -IMU2. data (: ,4) ; 

IMU_time = IMU2. data (: , 2) ; 

% Plot to show sectioned data 

plot (RCIN. data (; ,2) h: 10^-6,RCIN. data (: ,8) -1500,IMU2. data (: ,2) mUO^- 6,p 
:i:(180/pi) 

RCIN. data (: ,2) H:10^-6,RCIN. data (: ,9) -1500,IMU2. data (: ,2) MU0^-6,q 
Mc(180/pi)) 

legend ( ’ch6 ’ , 'p ’ , ’ch7 ’ , ’q ’ ) 
xlabel(’time (seconds)’) 
y 1 a b e 1 ( ’ pwm counts ’ ) 

[start_index , stop_index] = find_trigger (roll_trigger_data , 1700); 

f p ri n t f (’ Which test would you like to analyze? [Pick a number between 1 
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and %i]\n’ ,length(start_index)); 
num_of_test_desired = input('‘); 

%[ input, output, time ] = section_triggered_data ( input_data , 
input_data_time , output_data , output_data_time , trigger_data , 
trigger_time , num_of_test_desired ) 

[ input, output, time ] = section_triggered_data ( roll_command_data , 

command_data_time , p, IMU_time , roll_trigger_data , trigger_data_time 
, num_of_test_desired ) ; 

%Center controls for SysID and plotting 
input = input -1500; %center stick ouputs 1500 

% Theoretical Chirp 

fs = 1/50; 

t = [0;fs :7]; 

f0= 0.01;Wz 

fl= 10; Wz 

k = (fl-f0)/(t(end)); 

phi_0 = 0; 

input_theoretical = sin ( phi_0+2ii= pi fO^(7-t)+(k/2) - t) .^2)) ; 

input_theoretical = input_theoretical ^75; 

figure 

pi ot (time , input o ’ , time, output if:( 1 80/ pi ), t , i np u t_th e ore ti c al) 
legend(’roll cmd ’ , ’roll rate’) 

System Identification 
loop_rate = mean ( diff (IMU_time)) 10^-6; 
time_clean = [0: loop_rate :( length (time )-l)ii: loop_rate ]’; 
output_clean = interp 1 (time , output , time_clean pchip ’) ; 

%input_clean = interp 1 (time , input , time_clean ,’pchip ’) ; 
input_clean = interpl(t,input_theoretical,time_clean, ’pchip’); 


figure 

plot(time, input_clean , time, output_clean if:(l 80/pi )) 
legend(’roll cmd ’ , 'roll rate’) 
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[y,t,x] = estimate_model (input_clean , output_clean , time_clean); 
figure 

subplot (2,1 ,1) 

plot(time_clean ,input_clean) 

1 e g e n d ( ’ command ’ ) 
xlabel(’time (seconds)’) 
ylabel(’\delta (pwm) ’ ) 
title(’Roll Analysis’) 
subplot (2,1 ,2) 

plot (t ,y, time_clean , output_clean) 
legend ( ’ model ’ , ’ measured response ’) 
xlabel ( ’Time (s)’ ) 
ylabel(’Roll rate (rad/s)’) 
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APPENDIX G: 

Xi Adaptive Controller Source Code 


G.l £i Adaptive Control Source Code 

I* 

This program is free software: you can redistribute it and/or modify 
it under the terms of the GNU General Public License as published by 
the Free Software Foundation, either version 3 of the License, or 
(at your option) any later version. 

This program is distributed in the hope that it will be useful , 
but WITHOUT ANY WARRANTY; without even the implied warranty of 
MERCHANTABILITY or FITNESS EOR A PARTICULAR PURPOSE. See the 
GNU General Public License for more details. 

You should have received a copy of the GNU General Public License 
along with this program. If not, see <http://www.gnu.org/licenses/>. 

*/ 

I* 

Adaptive controller by Ryan Beall 
See 

http :/ / naira -hovakimyan . mechse . illinois .edu/ll-adaptive-control- 
tutorials / 

for mathematical basis 
*1 

#include <AP_HAL/AP_HAL. h> 

#include <GCS_MAVLink/GCS. h> 

#include " ADAP_Control. h" 

#include <stdio.h> 

extern const AP_HAL::HAL& ha 1 ; 

const AP_Param : : Groupinfo ADAP_Control :: var_info [ ] = { 

// adaptive control parameters 

AP_GROUPINPO_FLAGS("CH" , 1, ADAP_Control , enable_chan , 0, 
AP_PARAM_FLAG_ENABLE) , 

AP_GROUPINPO("AL" , 2, ADAP_Control , alpha, 24), 


89 





31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 
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62 

63 

64 

65 

66 

67 

68 

69 


AP_GROUPINFO( 

"GAMT" , 

3 , 

ADAP_ 

_Control , 

gamma_theta , 1000) , 

AP_GROUPINFO( 

"Q\MW" , 

4, 

ADAP_ 

_Control , 

gatnma_otnega, 1000), 

AP_GROUPINFO( 

"GAMS" , 

5 , 

ADAP_ 

_Control , 

gamtna_sigma, 1000), 

AP_GROUPINFO( 

"THU" , 

6, 

ADAP_ 

_Control , 

theta_max , 2.0), 

AP_GROUPINFO( 

"THL" , 

7, 

ADAP_ 

_Control , 

theta_min , 0.5), 

AP_GROUPINFO( 

"THE" , 

8, 

ADAP_ 

_Control , 

theta_epsilon , 5925) 

AP_GROUPINFO( 

"CMU" , 

9, 

ADAP_ 

_Control , 

otnega_tnax , 2.0) , 

AP_GROUPINFO( 

"OML" , 

10, 

ADAP_ 

_Control , 

omega_min , 0.5), 

AP_GROUPINFO( 

"OME" , 

11, 

ADAP_ 

_Control , 

omega_epsilon , 5925) 

AP_GROUPINFO( 

"SIGU" , 

12, 

ADAP_ 

_Control , 

sigma_max , 0.1), 

AP_GROUPINFO( 

"SIGL" , 

13, 

ADAP_ 

_Control , 

sigma_min , -0.1) , 

AP_GROUPINFO( 

"SIGE" , 

14, 

ADAP_ 

_Control , 

sigma_epsilon , 203) , 

AP_GROUPINFO( 

"WO" , 

15, 

ADAP_ 

_Control , 

wO, 25), 

AP_GROUPINFO( 

"K" , 

16, 

ADAP_ 

_Control , 

k, 0.45), 

AP_GROUPINFO( 

"KG" , 

17, 

ADAP_ 

_Control , 

kg, 1.0), 

AP_GROUPEND 







}; 

/* 

return true when enabled 
*/ 

bool ADAP_Control :: enabled ( void ) const 

{ 

// Enables Adaptive Controller instead of PID if rc PWM value is 
above 1700 milli seconds 

return (enable_chan > 0 && hal . rein->read (enable_chan-1) >= 1700); 

} 


/* 

reset to startup values 

*/ 

void ADAP_Control :: re s e t ( uint 1 6_t loop_rate_hz ) 

{ 

x_m = X; 
u = 0.0; 

u_lowpass = 0.0; 
u_sp = 0.0; 
theta = 1.0; 
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omega = 1.0; 
sigma = 0.0; 
integrator = 0.0; 

float u_cutoff_hz = wO/(2iJ:M_PI) ; //convert cutoff freq from rad/s to 
hz 

u_filter . set_cutoff_frequency(loop_rate_hz , u_cutoff_hz) ; 
u_filter . reset 0 ; 

float r_cutoff_hz = (alpha-2)/(2ii=M_PI) ; //convert cutoff freq from 
rad/s to hz 

r_filter . set_cutoff_frequency (loop_rate_hz , r_cutoff_hz ) ; 
r_filter . reset () ; 


1 ^ 

trapezoidal integration helper function 

float ADAP_Control :: trapezoidal_integration ( float yO , float yl_dot , 
float dt , float &y0_dot) 

{ 

float yl = yO + (dt/2)*(y0_dot + yl_dot) ; 
y0_dot = yl_dot; 

return yl; 

} 

I* 

adaptive control update. Given a target rate in radians/second and a 
current sensor rate on the same axis in radians / second , return an 
actuator value from -1 to 1 
*/ 

float ADAP_Control :: update ( uint 1 6_t loop_rate_hz , float target_rate , 
float sensor_rate , float scaler , float aspeed) 

{ 

float dt ; 

const float u_limit = radians (45); 

X = sensor_rate; 
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r = target_rate; 
r = r_fi1ter.apply(r); 

//reset reference model at initialization 
uint64_t now = AP_HAL: : micros64 () ; 

if (last_run_us == 0 || now - last_run_us > 200000UL) { 
reset ( loop_rate_hz ) ; 
last_run_us = now; 
return 0; 


dt = (now - last_run_us) ^ l.Oe-6; 
last_run_us = now; 


// u (controller output to plant) 

eta = thetaif:x + omegaH^u_lowpass + sigma; 

u_sp = eta; // u to state predictor 

u = constrain_float(eta-(kgii=r),-radians(90)/dt, radians(90)/dt); 

// Check Controller Saturation from previous time step 
bool saturated = ((u < 0 && u_lowpass >= 0.99^ u_limit) || 

(u > 0 && u_lowpass <= -0.99^ u_limit)) ; 

// kD(s) (cascaded second order low pass + simple integrator) 
u_lowpass = u_filter . apply (u) ; 

if (! saturated) { 

integrator = trapezoidal_integration( integrator , u_lowpass , dt , 
outl); 

u_lowpass = constrain_float (-k^i^integrator u_limit , u_limit) ; 

} 


// State Predictor (first order single pole recursive filter) 
// Reference/Companion Model 

float alpha_filt = exp(-alpha * dt) ; //alpha in rad/s 
alpha_filt = constrain_float ( alpha_filt , 0.0, 1.0); 

float beta_filt = 1-alpha_filt ; 
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159 

160 
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x_m 


= alph a_f i 11 ii=x_m + b e t a_f i 11 u_sp) ; 
x_error = x_m - x; 


// Constrain error to +-300 deg/s 

x_error = constrain_f1oat(x_error ,-radians(300), radians(300)); 


// Calculate solution to Lyapunov 1x1 matrix 
float Pb = l/(2:i: alpha) ; 


// Projection Operator 

theta_dot = proJection_operator( theta ,-gamma_theta:i:x_error^i^PbH^x , 
theta_epsilon , theta_max , theta_min) ; 
omega_dot = projection_operator (omega,-gamma_omegaH^x_error^f^PbH^ 
u_lowpass , omega_epsilon , omega_max , omega_min) ; 
sigma_dot = proj ection_oper ator ( sigma ,-gamma_sigmaH^ x_error ^f^Pb , 
sigma_epsilon , sigma_max , sigma_min) ; 

// Parameter Update using Trapezoidal integration 
if (! saturated) { 

theta = trapezoidal_integration (theta , theta_dot , dt , thetal); 

omega = trapezoidal_integration (omega, omega_dot , dt , omegal); 

sigma = trapezoidal_integration ( sigma , sigma_dot , dt , sigmal ) ; 

} 


theta = constrain_float (theta , theta_min , theta_max); 

omega = constrain_float (omega , omega_min , omega_max) ; 

sigma = constrain_float ( sigma , sigma_min , sigma_max); 

// Log Data to flash 

DataFlash_Class :: instance () ->Log_Write (log_msg_name , "TimeUS , Dt, 

A theta , Aomega , A sigma , Aeta , Axm, Ax, Ar , Axerr , AuL" , "Qffffffffff", 

now, 
dt , 

theta , 
omega, 
sigma , 
eta , 

degrees (x_m) , 
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degrees(x), 
degrees(r), 
degrees(x_error), 
degrees(u_lowpass)) ; 


return constrain_float ( u_lowpass / u_limit , -1, 1); 

} 


float ADAP_Control :: proj e cti on_oper ator ( f 1 o a t Theta, float y, float 
epsilon , float th_max , float th_min) const 

{ 


// Calculate convex function 

// Nominal un-saturated value is above zero line on a parabolic 
curve 

// Steepness of curve is set by epsilon 

float f_diff2 = (th_max-th_min) *( th_max-th_min ) ; 

float f_theta = (-4*(th_min - Theta) * (th_max - Theta))/( 
epsilon*f_diff2 ) ; 

float f_theta_dot = (4*(th_min + th_max - (2* Theta)))/( epsilon * 
f_diff2 ) ; 

float proj ection_out = y; 

if (f_theta <= 0 && f_theta_dot*y < 0) 

{ 

projection_out = y*(f_theta+1); //y-(y*(-1 *f_theta)); 

} 


return 

} 


projection_out ; 


/* 

send ADAP_TUNING message 

*/ 

void ADAP_Control ;; adaptive_tuning_send ( mavlink_channel_t chan, uint8_t 
axis ) 


if (! enabled!) || !HAVE_PAYLOAD_SPACE( chan , ADAP_TUNING)) { 
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return ; 


} 

mavlink_msg_adap_tuning_send ( chan , axis , 

r , 

X, 

x_error , 
theta , 
omega, 
sigma , 
theta_dot , 
omega_dot , 
sigma_dot , 
x_m, 

u_lowpass , 

u); 
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